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b rib width 

c c £l , c a turbulence model constants 

c, coefficient in source term of dissipation rate equation 

c p constant pressure specific heat 
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d x inner diameter 
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d e d 2 — d ] 

dj jet diameter 

E empirical constant in the law of the wall 

e roughness height 

/ friction factor 
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e + roughness Reynolds number, — — 

k turbulence kinetic energy, k = u ' u ! 
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mass flow rate 
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pressure 
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volumetric flow rate 
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radial coordinate 
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inner radius of the annulus 
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outer radius of the annulus 

R 

radius of pipe 

4 tu 

Re D 
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Re„, 


w 

Pw Pb 

Si 

source term for <p t 
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T 

local temperature 
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dimensionless temperature, pu c p (T w - T)j q w 

u 
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u* friction velocity, 

u + dimensionless velocity, 

u 

u~uj Reynolds stress 

y+ dimensionless distance from wall, 

w rib width 



yu_ 
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Greek symbols 

a, p constants for algebraic model 

£ turbulence dissipation rate 

£ modified dissipation variable 

T thermal conductivity 

p dynamic viscosity 

p t turbulent viscosity 

v kinematic viscosity 

v, turbulent kinematic viscosity 

a £ turbulent Prandtl number for </>,</> = k, e, h 

p density 

c fry pressure strain tensor 

V volume 

ip level of turbulence intensity 

Subscripts 

av average 

b bulk 


vm 



c constant property 

in inlet 

P grid node at center of control volume 

t turbulent 

v viscous sublayer 

w wall 

e,w,n,s control volume faces 

£, grid nodes at east, west, north, and south of control volume 
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ABSTRACT OF THE DISSERTATION 


COMPUTATION OF TURBULENT RECIRCULATING FLOW 
IN CHANNELS, AND FOR IMPINGEMENT COOLING 

by 

Byong Hoon Chang 

doctor of philosophy in Mechanical Engineering 
University of California, Los Angeles, 1992. 

Professor Anthony F. Mills, Chair 

Fully elliptic forms of the transport equations have been solved numerically 
for two flow configurations. The first is turbulent flow' in a channel with 
transverse rectangular ribs, and the second is impingement cooling of a plane 
surface. Both flows are relevant to proposed designs for active cooling of 
hypersonic vehicles using supercritical hydrogen as the coolant. Flow down- 
stream of an abrupt pipe expansion and of a backward-facing step w'erc also 
solved with various near-w'all turbulence models as benchmark problems. A 
simple form of periodicity boundary condition w 7 as used for the channel flow 
with transverse rectangular ribs. The effects of various parameters on heat 
transfer in channel flow with transverse ribs and in impingement cooling w'ere 
investigated using the Yap modified Jones and Launder low Reynolds number 
k - e turbulence model. For the channel flow, predictions were in adequate 
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agreement with experiment for constant property flow, with the results tor 
friction superior to those for heat transfer. For impingement cooling, the 
agreement with experiment was generally good, but the results suggest that 
improved modelling of the dissipation rate of turbulence kinetic energy is re- 
quired in order to obtain improved heat transfer prediction, especially near the 
stagnation point. The k — e turbulence model was used to predict the mean 
flow and heat transfer for constant and variable property flows. The effect of 
variable properties for channel flow was investigated using the same turbu- 
lence model, but comparison with experiment yielded no clear conclusions. 
Also, the wall function method was modified for use in the variable properties 
flow with a non-adiabatic surface, and an empirical model is suggested to cor- 
rectly account for the behavior of the viscous sublayer with heating. 
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Chapter I 

INTRODUCTION 


1.1 BACKGROUND 

Current interest in hypersonic flight has led to renewed activity related to 
high temperature structures. A critical problem is to ensure the survivability 
of components subjected to intense aerodynamic heating, such as the nose, 
wing leading edges, and the engine inlet. Hydrogen-fueled scramjet engines are 
under development, and a cooling system which shows considerable promise 
is based on the use of the hydrogen fuel as a coolant before it is injected into 
the combustor. Engine inlets require panels in which the hydrogen flows 
through channels beneath the skin. The use of enhanced surfaces to increase 
heat transfer coefficients inside the channels is an attractive option to improve 
performance. For the nose, or wing leading edges, impingement cooling is a 
possible approach. A jet of hydrogen impinges on the backface of the skin at 
the stagnation point or line, and flows rearwards. 

Experimental data will certainly be required in order to develop systems 
described above. But test work using supercritical hydrogen is both expensive 
and hazardous: thus the use of modern computational fluid dynamics(CFD) 
methods are an attractive partial alternative. If successful computer models 
can be developed, the scope of the experimental program can be reduced ac- 
cordingly. The CFD model then becomes a tool for interpolation in, and 
modest extrapolation of, the experimental data. 
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1.2 OBJECTIVES OF THE PRESENT STUDY 

There are numerous problems relevant to active cooling, and suitable for 
CFD modeling. Two problems have been chosen for the present study. These 
arc : 

1. The development of roughness functions for transverse ribs used to 
augment heat transfer in channel flow. These functions will facilitate 
the engineering calculation of pressure drop and heat transfer in cooling 
panels. 

2. Prediction of the flow field and heat transfer rates for two- dimensional 
impingement cooling. The results will facilitate the design of the 
impingement system as a heat exchanger. 

These problems have three important features in common. Firstly, both 
involve recirculating flows with flow separations and reattachments. Thus the 
elliptic form of the transport equations must be solved. Secondly, both flows 
are turbulent, and due to the critical importance of heat transfer near flow re- 
attachment points, simple turbulence models, such as the mixing length model, 
are inadequate. Thirdly, fluid property variations are large adjacent to the 
wall of primary interest. It is therefore necessary to solve the conservation 
equations in the near wall region: so called conventional wall functions should 
not be used to bridge the region. Hence an essential requirement of this study 
is a suitable turbulence model which can be used for recirculating flows with 
large property variations. 

A mixing-length turbulence model is not suitable for these flows because the 
turbulent viscosity and thermal conductivity vanish where the mean velocity 
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gradient is zero; also, the influence of convection and diffusion on the turbu- 
lence kinetic energy are not accounted for. Two-equation models of turbu- 
lence arc more suitable for these flows, and many variations have been used 
by prior workers. The conventional wall function approach replaces the near 
wall region with formulas based on the logarithmic velocity and temperature 
profiles in order to avoid using a fine grid near the wall. However, unless 
special wall functions are developed, which properly account for variable 
property effects, a fine grid is required near the wall. 

Since experimental data for the problems of interest is sparse, related prob- 
lems, namely heat transfer downstream of a sudden pipe expansion and of a 
backward-facing step, will be used as benchmarks to initially test the turbu- 
lence model, and the computational technique. These flows have been the 
subject of many experimental and numerical studies: there is both a good data 
base, and a well documented history of turbulence model development. 


1.3 LITERATURE SURVEY 

1.3.1 Flow downstream of a sudden pipe expansion 

The flow over a backward-facing step or the flow downstream of a sudden 
pipe expansion is very complex. It involves recirculating flow, shear layer re- 
attachment, a counterrotating secondary vortex in the corner, and boundary 
layer redevelopment. According to Eaton and Johnston [1], the length of the 
separation region behind the step fluctuates so that the impingement point of 
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the separated shear layer is not stationary. The measured maximum backward 
velocity in the recirculating flow region is reported to be usually over 20 % of 
the freestream velocity. The k — e two-equation model that determines the 
turbulent viscosity from the solution of two transport equations, namely the 
turbulent kinetic energy, k, equation and its dissipation rate, c, equation, has 
been used by various workers to solve this recirculating flow. The standard 
model is applicable only to regions of high Reynolds number, and its use with 
wall function approach to bridge the viscous sublayer has been popular in or- 
der to reduce computing costs. But, a conventional wall function based on the 
friction velocity, , is not appropriate for this flow since the predicted 

heat transfer is zero at separation and reattachment points. This deficiency in 
the wall function can be removed by adopting the turbulence kinetic energy 
as a velocity scale, as proposed by Launder and Spalding [2]. A wall function 
approach based on constant non-dimensional viscous sublayer thickness was 
attempted by Chieng and Launder [3]. A coding error of Chieng and 
Launder was detected by Johnson and Launder [ 4], and a subsequent recal- 
culation showed underprediction of Nusselt number for the experimental heat 
transfer data of Zemanick and Dougall [ 5]. Johnson and Launder [4] ob- 
tained good results by making the nondimensional viscous sublayer thickness 
a linear function of the ratio of the rate of turbulence energy diffusion to the 
rate of turbulence energy dissipation in the sublayer. Amano [ 6] extended 
this two-layer model to a three-layer model by introducing a buffer layer be- 
tween the viscous and turbulent regions. The prediction by the three-layer 
model compared favorably with the experimental data of Amano et al. [7], 
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Launder [ 8] has made a critical evaluation of the wall function approach, and 
concluded that a more appropriate form of the turbulence energy dissipation 
rate equation needs to be established before further refinement of wall func- 
tions. Also, the use of wall functions failed to predict the secondary corner 
vortex. Baughn ct al. [9] have made extensive experimental measurements 
of the local heat transfer coefficients to an air flow downstream of an 
axisymmctric abrupt expansion in a circular pipe with a constant wall heat 
flux. The runs were made with the expansion ratio from 0.267 to 0.8 and over 
the Reynolds number range of 5300 to 87000. The maximum Nusselt number 
was almost eleven times larger than that for fully developed pipe flow, as 
given by the Dittus-Boelter relation, for an expansion ratio of 0.266 based on 
pipe diameters, and a downstream Reynolds number of 8,112. Launder [10] 
presented a review on the various methods for computing heat transfer coeffi- 
cients in complex turbulent flows. A significant improvement in the calcu- 
lation of heat transfer rates by Yap [11] was the addition of a source term to 
the dissipation rate equation to reduce the excessive turbulence near-wall 
length scale. A comparison with Baughn's experimental data showed that the 
method of adopting a low Reynolds number k — £ model across the sublayer 
and an algebraic stress model beyond gave the best heat transfer prediction 
with a good Reynolds number dependence. However, the prediction by a low 
Reynolds number k - e model with the Yap's correction was fairly good. With 
the low Reynolds number k-e model, computations are performed all the 
way to the wall, and the algebraic stress model employs a simplified form of 
the Reynolds stress transport equations. 
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1.3.2 Flow downstream of a backward-facing step 

Aeronautical researchers have long been interested in the flow over a 
backward-facing step because of the practical importance of predicting base 
pressure of bluff bodies moving with high speed(such as bullets and coasting 
missiles) [12]. More work has been done in supersonic flow, but the low speed 
flow over a backward-facing step is also used as a building block flow for 
workers developing turbulence models [1]. Some of more relevant work in 
literature to the present study are listed here. 

Eaton and Johnston [1] were among the first to make measurements in the 
highly unsteady and reversing flow behind a backward-facing step using a 
pulsed-wire anemometer, thermal tuft, and pulsed-wall probe. Turbulence 
quantities and the skin friction were measured, and it was found that the tur- 
bulence intensity decreases rapidly downstream of reattachmcnt in the rede- 
veloping boundary layer with the turbulence actually beginning to decay 
upstream of reattachment some one or two step heights. Durst and Tropea 
[13] used a water channel to study a backward-facing step flow, acquiring the 
data with a laser-doppler anemometer. In their study, they varied the expan- 
sion ratio and found a strong dependence of the rcattachment length on this 
quantity. Driver and Seegmillcr [14] acquired wind-tunnel fluid dynamic data 
for the backward-facing step using a laser-doppler anemometer. They ob- 
tained many of the important turbulence quantities, such as production and 
dissipation rates in the flow. In their comparison with various numerical 
methods, they found that a modified algebraic stress model predicts the flow 
more accurately than does the k - £ model or the unmodified algebraic model. 
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The modified algebraic stress model was found to be more sensitive to 
streamline curvature effects on the turbulence. 

Scban [15] measured velocity and temperature distribution downstream of 
a backward facing step using a constant heat flux surface. For a one-inch step 
and a free-stream velocity of 150 ftf sec, he found that the mean reattachment 
point was about six step heights downstream of the step. The typical law of 
the wall was not found in either the separation or rcattachment regions. Filctti 
and Kays [16] measured the local heat transfer coefficient on a constapt- 
temperature surface behind a symmetric sudden expansion. The expansion 
ration of 2:1 was used, and the peak heat transfer rate occured at the reat- 
tachment point, which was approximately four step heights downstream of the 
step. Significant augmentation of heat transfer rate downstream of reattach- 
ment, up to six times the fiat-plate value, was found. However the measure- 
ments were taken from the step edge to 14 step heights downstream and the 
Nusselt number profile did not approach the typical flat-plate value. The peak 
heat transfer rate based on Nusselt number was found to vary as Re 0 6 over the 
range 70,000 < Re < 205,000. Seki et al. [17], [18] obtained temperature and 
velocity measurements in the flow over the constant heat flux surface of a 
double-sided, backward-facing step. The velocity and temperature measure- 
ments were correlated to measure v'T , the turbulent transport of energy in the 
direction normal to the wall. The v'T* profiles indicated that the turbulent 
transport of heat in the direction normal to the surface increases as Re 2/3 and 
is proportional to the heat flux rate from the wall. Vogel and Eaton [19] used 
a single sided sudden expansion with a constant heat flux surface to make 
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combined heat transfer and fluid dynamic measurements in a separated and 
reattaching boundary layer. Stanton number profiles were obtained for four 
different Reynolds numbers ranging from 13,000 to 42,000, and the bottom 
wall of the development section was porous to vary the boundary layer thick- 
ness at the test section entrance from 0.3 to 6 cm. The upstream boundary 
layer thickness had a significant effect of the heat transfer rate near reattach- 
ment but very little effect either up- or downstream of reattachment. The 
temperature profiles showed that the heat transfer resistance is dominated by 
the near-wall region. 

Sindir [20] performed a numerical study of the effects of expansion ratio 
on two-dimensional separating and reattaching flow in plane backward-facing 
step with four models of turbulence. The k — e model, modified k — e model, 
algebraic stress model, and modified algebraic stress model were used, and the 
modified versions employed a production term in the dissipation equation that 
was made more sensitive to streamline curvature effects. The modified alge- 
braic stress model produced the best predictions in the reverse flow region but 
performed more poorly than others in the redevelopment region. Heat transfer 
calculations were not performed. Gooray et al. [21] employed a two-pass 
procedure: the first pass with an improved k — e model with the standard wall 
function to find the reattachment point and the second pass with an improved 
low- Reynolds number model applied to downstream of the reattachment. 
Good prediction of the local Nusselt number was obtained downstream of the 
step, but comparasion was made with limited data. Scherer and Wittig [22] 
have used the k — e model with one-layer and two-layer wall function ap- 


8 



proach of Chicng and Launder [3]. The rcattachment length was undcrpre- 
dicted by 10 to 20% for the Reynolds number from 40,000 to 85,000. The 
peak heat transfer was significantly undcrpredictcd by the one-layer model, 
and the two-layer model gave good prediction of local Nusselt number down- 
stream of the step. Ciofalo and Collins [23] have noted that the method of 
Johnson and Launder [4] involves the near-wall profile of the turbulence 
kinetic energy, and is sensitive to numerical errors. They have noted that 
slopes of the nondimensional velocity profiles in the reverse flow region and 
shortly downstream of reattachment point were not far from the equilibrium 
value, while the intersection with the linear region was greatly reduced. They 
modified the standard wall functions to relate the non-dimensional viscous 
sublayer thickness to the near-wall turbulence intensity. Comparison of the 
Nusselt number prediction with the experiment data on backward-facing step 
showed significant improvement over the standard wall functions. However, 
the prediction downstream of the reatttachment point converged to a value 
about 12% lower than that of the experimental data. 


1.3.3 Roughness functions for transverse ribs 

The pioneering work on roughness functions was done by Nikuradse [ 24] 
who performed extensive experimental measurements on circular pipes with 
sand of a definite grain size glued on the inside surface. He showed that the 
dimensionless velocity distribution is given by : 

u + = 2.5 In (-y) + R (e + ) (1) 
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where R (e + ) is termed the roughness function. 

A similar approach based on a modified Reynolds analogy was used by 
Dipprey and Sabersky [25] to correlate their heat transfer data for flow inside 
tubes with a sand grain indentation surface. The Stanton number was given 

as : 


c ± 

2 

\+J^lH(e + ,Pr)-R(e + )l 


( 2 ) 


where the heat transfer roughness function H ( e + , Pr) was correlated as 


, °- 22 o 44 
// = 5.19 k+ Pr 


( 3 ) 


Following Nikuradse, R (e+) was taken as 8.5 for fully rough flow. 

Webb, Eckert, and Goldstein [26] used their experimental data to develop 

the roughness functions R and H for flow in a tube of repeated transverse rib 

roughness. For 0.01 < -j- < 0.04, 10 < < 40, and e + > 35, 

D t 


, 0.53 

R =0.95(^-) 


(4) 


0.28 

H = 4.5 e + 


Pr 


0.57 


(5) 


Han et al. [27] performed experiments to study additional effects of rib 
shape and angle of attack of ribs to main flow in a channel with repeated 
transverse rib roughness. A 45“ flow attack angle was found to give higher 
heat transfer for the same pressure drop, as compared to a 90“ flow attack 
angle. Their roughness functions for 90“ flow attack angle are 
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( 6 ) 


R = 0.97 (-J-) 0 53 for LI e> 10, e + > 35 
= 4.45 ( A )" 013 for L I e < 10, e + < 35 

H = 5.05 e + ° 2 Vr 0 - 57 for e + > 35 (7) 

Han ct al. based the heat transfer coefficient on the total area of the heat 
transfer surface including the rib area, for the ribs constituted an appreciable 
fraction of the total area for a small value of L/e. When the heat transfer co- 
efficient is based on the projected area, the heat transfer roughness function 
becomes 

H = 4.04 £> + ° 28 />r° 57 for e + > 35 (8) 

Dalle Donne and Meyer [28] performed experiments in an annulus with an 
inner surface of transverse rectangular ribs and a smooth outer surface. For, 
e + > 30, the roughness functions were correlated as 

R = R (oo) + 0.4 In ( e — — ) (9) 

V ' 0.01D/2' 

where 

fl(oo) = 9.3(-^-r 073 - [ 2 + (L _ 7 fc)/(| ] log | 0 (y) (10) 

for 1 < L ~ — < 6.3 
e 

R( oo) = 1 .04 (^=A) 0 - 46 - [2 + (L _ 7 6)/e J lQ glo(f) 01) 

for 6.3 < L ~- < 160 
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H = (4.16 e 


e 


+0.282 


+ 


53 

TIT 


■)/>r°' 44 (2+) 05 [- 


0.01 (r 2 -r,) 


] 


0.053 


( 12 ) 


The first attempt to analytically determine the roughness functions for flow 
over rectangular ribs was made by Lewis [29]. The flow was approximated 
by a series of attached and separated flow regions, and some empirical infor- 
mation from experiments over cavities and steps was required. The A: - s tur- 
bulence model with the wall function boundary condition was used by Lee et 
al. [ 30] to predict roughness functions in an annulus with ring type rectan- 
gular roughness on the inner pipe. In a numerical study, fully developed flow 
in a single module was solved using the periodicity conditions, as proposed by 
Patankar et al. [31], in order to avoid the entrance region problem. A cor- 
rection to the turbulent viscosity similar to that of Leschziner and Rodi [32] 
was used to account for extra strain rates due to streamline curvature. How- 
ever, as shown by Launder [10 ] the use of a wall function gives rather poor 
prediction of Nusselt number for flow in an abrupt pipe expansion and flow 
around a 180° square-sectioned bend. 


1.3.4 Variable properties flow 

There have been numerous experimental studies on heat transfer in turbu- 
lent flow in smooth round tubes with large wall to bulk temperature ratio, and 
many correlation methods have been proposed. A good review is given by 
Petukhov [33]. For heat transfer, a recommended correlation is 
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where C and n arc constants far from the entrance. Petukhov observes that 
different investigators obtained different values for n depending on the range 
of TJ T b . For variable property gases, he suggests 

n = -{a log (-=r~) + 0.36) (14) 

J b 


a = 0 for cooling, and a = 0.3 for heating was recommended. 


For the friction factor, the recommended correlation is 




m 

) 


(15) 


where 

m = —0.6 + 5.6 Re w ®'^ (16) 

Petukhov used the Reichardt eddy diffusivity profile to obtain analytical re- 
sults. 

Sleicher and Rouse [34] found a better fit to large amount of heating data 
when n was modified to 

T 1/4 

n = — log (— ~~) +0.3, 1 < TJT b <5, x/Z) > 40 (17) 

1 b 

Many correlations are of similar form and three that are among the many 
reviewed by Petukhov are summarized below. McEligot et al. [35] have per- 



formed experiments with air, nitrogen, and helium in a tube with entering Re 
from 15,000 to 233,000. The following forms were recommended 


0.4 T„ -0 5 


Nu b = 0.021 ReFPrF(-f-) 


C I + (ff 


- 0 . 7 - 


w 


fb_ _ , 

fc T b 


T - 0.1 


(18) 


) 


for I < TjT b < 2.5 


Lelchuk and Dyadyakin [36] recommend the same correlation for the heat 
transfer without the entry length correction term, but for the friction factor the 

exponent on the temperature ratio, m, was recommended as 

T 

- 0.16 for 1.3 < — 1 — < 2.4. Perkins and Worsoe-Schmidt [37] used precooled 

Tb 

nitrogen to obtain local values of TjT b from 1.24 to 7.54. Recommended 
correlations were 


T - 0.7 - 0.7 T 0.7 

Nu b = Re b % Pr b A {-jT) [ 1 + (-^-) (■yp 3 

f T - 0.6 

J_w_ _ , _w , 

Jc~ T b 


(19) 


Perkins found that when the correlation of friction factor was tried on the bulk 
Reynolds number, the exponent, n, had to increase from 0.1 to 0.3 over the 
temperature ratio range of the experiment. The above correlation for friction 
factor was based on the wall Reynolds number. 

Experiments with repeated-rib roughness were performed by Vilemas and 
Simonis [38]. Air was used in annuli with a rough inner wall for 
5 x 10 3 < Re < 5 x 1 0 5 , 1 <TjT in < 2.8, 0.0028 < e\d e < 0.021, and 8.3 
<Lje< 13. For rectangular roughness djd 2 was 0.42, and for rounded 
trapezoidal roughness djd 2 was 0.35. The outer wall was smooth, and only the 
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inner wall heated. It was found that, unlike smooth channels, the influence 
of TJT b on the local heat transfer in rough annular channels depends heavily 
on Re. The effect of TJT b decreased with increasing Re. The influence of 
Re on the exponent, n, was greater for larger ejd e . For example, for a channel 
with ejd e = 0.021, n changed from about - 0.36 at Re = 2.8 x 10 4 to - 0.1 at 
Re = 3.5 x 10 5 . The exponent, n, was correlated within ± 10 % by 

„ = -(0.29 + 0Me il '^)Re 2 ‘ >el ‘ , '0 (20) 

The friction factor for annular channels with rectangular roughness was cor- 
related within ± 4 % by 

f b = (0.053 + 1.83 -j-)Re b Qm (21) 

In their experiment, an influence of variable properties on the friction was not 
observed in smooth channels for TjT in up to 1.8. For the rough channel, the 
effect of variable properties decreased with an increasing Re and became neg- 
ligible at Re of about 5 x 10 s . 

Wasscl and Mills [39] performed numerical calculations for variable prop- 
erties turbulent flow with cooling for both smooth and rough walls. For flow 
in rough pipes, the Nikuradse mixing length expression was used with a 
roughness form drag coefficient and sub-layer Stanton number characterizing 
transport to the wall. The variable properties in the roughness functions were 
evaluated at the characteristic roughness height. Flow in a pipe with the 
sandgrain roughness size equal to kJR = 15 and 60, and rectangular rib 
roughness with ejR = 0.02 and Lje = 10 were investigated. The results 
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showed that f b jf c depends on the roughness pattern and size, and a slight effect 
of Reynolds number was observed. For the large sandgrain size, the Reynolds 
number effect on St„ jSt c was large. The following correlations were recom- 
mended for the ribs. 


A A - 0 - 2 £!l - f^r 0 ' 25 
fc V ’ Sr c h b 


( 22 ) 


1.3.5 Jet impingement cooling 

There have been numerous experimental studies of both turbulent free jets 
and turbulent impinging jets in literature: some of the work relevant to the 
present study is described. Poreh, Tsuei, and Cermak [40] made measure- 
ments of mean velocities, turbulence intensities, Reynolds stresses, and the wall 
friction in a radial jet formed by an impinging circular jet on a smooth flat 
plate. The ranges of parameters were Hjd = 8 - 24, and Re — 64,000 - 
288,000, where H is the distance from jet outlet to the impingement surface 
and d is the jet outlet diameter. Beltaos and Rajaratnam [41] performed ex- 
periments on plane turbulent impinging jets over the range, Hjd = 14.04 - 67.5 
and Re = 5,270 - 9,400. Mean velocities, static pressure, and shear stresses 
over the impingement surface, and mean velocities and turbulent shear stresses 
over the free jet region have been reported. Beltaos and Rajaratnam [42] also 
made similar measurements for circular turbulent jet impingement over the 
range of Hjd = 21.2 - 65.7 and Re = 35,200 - 80,400. Giralt et al. [43] made 
measurements for circular turbulent impinging jets over the range, Hjd —1-2 
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- 25 and Re = 30,000 - 80,000. Velocity and length scales of the impingement 
flow field were used to scale impinging jet centerline velocities and pressure 
distributions. Wygnanski and Fiedler [44] made measurements of mean ve- 
locities and turbulence quantities for an axisymmetric turbulent free jet. It was 
concluded that the jet was truly self-preserving some 70 diameters downstream 
of the nozzle. The Reynolds number was in the order of 100,000, and most 
of the measurements were made for the distance of 40 to 100 diameters 
downstream. 

Gardon and Akfirat [45] were among the First workers to make thorough 
heat transfer measurements for two-dimensional turbulent impinging jets. 
Their results showed that stagnation heat transfer coefficients can be increased 
by artificially increasing the initial turbulence of the jet, with the effect being 
the largest for Hjd < 8. A secondary peak in heat transfer coefficient for 
Hid = 2 gradually disappeared with increasing induced turbulence. With 
circular jets, they also observed the peak heat transfer coefficient at r= 1/2 d 
(rather than at the stagnation point) over a range of nozzle-to-plate spacings, 
up to about Hjd = 3. Goldstein and Behbahani [46 ] performed experiments 
to study the effects of cross flow on a circular jet impinging on a flat plate. 
The maximum Nusselt number was found to decrease with increasing cross 
flow for jet-to-plate spacing, Ljdj , of 12, but for Ljdj = 6, the maximum 
Nusselt number increased with moderate cross flow when (p J u J )l(p oc u oc )> 9. 
Hrycak [47] made stagnation point heat transfer measurements for round jets 
impinging on a flat plate, and found a Nusselt number dependence on the 1/2 
power of the Reynolds number over the range, Hjd =1-20 and Re = 14,000 
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- 67,000. More recently, Baughn and Shimizu [48] made a careful measure- 
ment of heat transfer coefficients of a single air jet impinging on a flat surface. 
For Hjd = 2, a second maximum heat transfer occurred at rjd of approxi- 
mately 2. The local Nusselt number decreased to 70% of the Nusselt number 
at the stagnation point at about rjd = 1.5, increased to 80% at about rjd = 
2, and then decreased to 20% at rjd = 9. The experimental data for only one 
Reynolds number equal to 23,750 was reported, and the stagnation Nusselt 
number was 140 for Hjd = 2. As with other investigators, the maximiyn 
stagnation point heat transfer occurred at a Hjd of approximately 6. Unlike 
Gardon and Akfirat's measurements, the peak heat transfer rate occurred at 
the stagnation point, and not at r — 1/2 d for Hjd = 2. 

Wolfshtein [ 49] employed a one-equation turbulence model with wall 
function boundary condition for a jet normal to a flat plate. Only momentum 
equations were solved in the numerical study. Amano and Brandt [50] per- 
formed a numerical study of the flow characteristics of a turbulent jet with the 
high Reynolds number form of k - e model. Upper plate surface pressure, 
ground plane surface pressure, and velocity field for two-dimensional jet 
impingement were computed by Chuang [51] using the high Reynolds number 
k — e model. Similar computations were also performed by Hwang and Liu 
[52] with a fully developed jet velocity profiles along the inlet of computation 
domain. Agarwal and Bower [53] used the Jones and Launder's low- 
Reynolds number model to study the pressure distributions on the ground 
plane in the presence of upper surface in normal impingement of the 
compressible jet as in the case of VTOL aircraft. Good agreement was ob- 
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taincd. Malin [54] has used variations of the standard k — e model and 
k - W model of Spalding [55] to compute free jets and wall jets with the for- 
ward marching-integration procedure of PHOENICS. Inclusion of the 
irrotational strain terms in the production term of the turbulent kinetic energy 
generally gave improvements in the jet spreading rates. Predictions of the 
turbulent shear stress and the turbulent kinetic energy improved close to the 
center of the jet or the impingement wall, and somewhat deteriorated in the 
outside regions. 

A plane turbulent jet impinging obliquely at 70° on a flat surface was solved 
by Hwang and Tsou [56] using a two- equation turbulence model with wall 
function boundary conditions. The distribution of the inlet profile was ob- 
tained from the solution of plane turbulent jet, and the Nusselt number pred- 
ictions compared well with experimental data, except at the stagnation point 
where there seems to be about 32% overprediction. Rodi and Scheuerer [57] 
used Lam and Bremhorst's [58] low-Reynolds number version of the k - 1 
model to predict the heat transfer coefficients around gas turbine blades. The 
stagnation point was avoided by prescribing inlet profiles of dependent vari- 
ables. Predicted and measured heat transfer coefficients were in good agree- 
ment except in the transition region due to its short length predicted by the 
turbulence model. Computations of two-dimensional turbulent free jet and 
turbulent impinging jet were made by Looney and Walsh [59]. Generally 
good predictions of hydrodynamic and turbulence quantities were obtained 
with the standard k — £ model, but the results of the heat transfer coefficients 
on the impingement surface were generally poor. The algebraic stress model 
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gave improvement in prediction of both mean and turbulence quantities for a 
developing plane free jet. Amano and Sugiyama [60] employed Chieng and 
Launder's wall functions to predict heat transfer coefficients for an 
axisym metric jet impinging on a flat plate. A fairly good prediction for Hid 
= 4 and Re = 20,000 was obtained , but there was a 35% overprediction at 
the stagnation point, and about 27-33% underprediction over the range rid 
= 3-15 occured when compared with data for Hjd = 10. Polat et al. [61] 
studied the effect of various wall functions on prediction of heat transfer for a 
confined two- dimensional turbulent air jet impinging on a flat surface. Chieng 
and Launder's wall functions approach with the turbulent kinetic energy 
evaluated at the first node instead of at the edge of the viscous sublayer in the 
evaluation of wall shear stress gave the best prediction. Even though this 
method correctly predicted the secondary peak in heat transfer coefficient for 
Hjd = 2.6, the prediction became poor downstream of the stagnation point, 
with about 37% overprediction at y/d - 16. For Hjd > 6, the model contin- 
ued to predict off-stagnation minima and maxima even though such features 
are not present in experiment, thus the prediction became poor in the vicinity 
of the stagnation point also. As with Amano and Sugiyama, the wall function 
approach seems to perform poorly in the redeveloping region. Polat et al. 
[62] made flow and heat transfer predictions for confined turbulent impinging 
slot jets, with and without through-flow. Chieng and Launder's wall functions 
and the modified shear stress expression to account for the effect of mass 
transfer at the impingement surface were used. Generally good predictions for 
heat transfer coefficients within 10% were obtained for small through-flow, 
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but substantial overprediction at the stagnation point and at higher through- 
flow rates occured. Yap [1 1] has made calculations for the impinging jet flow 
experiment of Goldstein and Bchbahani [46], with both the algebraic stress 
model and the low Reynolds number k - e model. Comparisons of the two 
methods were made for predictions of centerline velocity decay and mean ve- 
locity and temperature profiles downstream of the stagnation point, and there 
was little difference. Local heat transfer predictions by both models were al- 
most the same, and both overpredicted the stagnation point heat transfer by 
15 %. Predictions were not in good agreement within about a nozzle diameter 
of the stagnation point, but improved further away. Only one value of the 
parameter H\d{= 6) and one value of Reynolds number was investigated. 

Hrycak [63] reports that there are a substantial number of investigations 
of heat transfer from jets impinging on flat plates but only relatively few ex- 
perimental results concerning heat transfer from jets impinging on concave 
surfaces. Some of the early experimental studies related to the concave inner 
surface of gas turbine airfoils cooled by impinging air jets were done by 
Metzger et al. [64 ], Chupp et al. [65], and Jusionis [66]. More relevant to 
the present study is the experimental work by Livingood and Gauntner [67] 
[68], [69 ]. Correlations of Nusselt number were presented in terms of the 
dimensionless quantities involving nozzle diameter, nozzle-to-target separation 
distance, target cylinder diameter, and nozzle center-to- center spacing for 
number of air jets more than one. The nozzle-to- target separation distance 
was found to have a greater effect on the normalized Nusselt number distrib- 
ution than did the variation of Reynolds number. Hrycak [63] performed a 
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similar experiment with a much more careful study of the hydrodynamics of 
the impinging jet. His heat transfer correlation was based on Froessling's 
hydrodynamics solution [70] and Colburn's analogy. 
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Chapter II 

MATHEMATICAL MODEL 

2.1 MEAN TRANSPORT EQUATIONS 

The set of elliptic partial differential equations governing a single- phase, 
compressible, variable property flow is as follows. 


Continuity'. 
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Momentum : 
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where the stress tensor, t„, is given by 


d^ duj 2 du k 


(26) 


23 



The instantaneous components in the above equations can be decomposed into 
sums of mean and fluctuating parts. For example, the streamwise velocity 
component, u, can be written as, 

u=U +u’ 

Ignoring fluctuations of physical properties except density, and time-averaging 
the decomposed equations yields the following. 

Continuity. 
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The equations (27), (28), and (29) cannot be solved in their exact form. Mod- 
elling of certain terms in the equations such as u-uj are required to provide a 
closed set of equations. 


2.2 TURBULENCE MODELS 

For steady flow and neglecting turbulent correlations involving density fluctu- 
ations, the transport equations reduce to the following. 

-j-(pUi)= 0 (30) 

CXj 


dx : 
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where the overbar notation for mean values and the prime notation for the 
fluctuating component of velocities has been dropped. The time mean of the 
product of the fluctuating velocities are modelled using the modified 
Boussinesq concept 
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The turbulent heat fluxes are approximated using the eddy diffusivity concept 


as 


— 77 f*-i dh 

y J Pr t dxj 


(34) 


2.2.1 The high- Reynolds number k-z model 

For the k — e model, the turbulent viscosity, p , , is modelled as 

(35) 


where c is a constant for high Reynolds number flow. 

The transport equation for turbulent kinetic energy k can be obtained from 
manipulation of the three normal stress equations [71 ]. 
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Following Jones and Launder [72], the diffusion term is modelled as 




d (Pi dk \ 
dxj ^ ° k dxj J 


(37) 


Substitution of the above equation into Eq.(36) and rewriting the last two 
terms on the RHS results in the following form valid for high Reynolds num- 
ber flow, 
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J-(pU j k) = ^-l-^ L -j L ) + P-pe 

dxj V J oxj \ °k dxj J 


where £ is the dissipation rate of turbulent kinetic energy defined as 


£ = V 


duj du i 
dxj dxj 


and P denotes the generation rate of turbulence energy. 


P=- P u i u j-— 

J h v . 


The transport equation for £ is obtained by taking the derivative of the in- 
stantaneous x t direction momentum equation with respect to Xj, multiplying 
by 2v du'jdxj and then averaging. 
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Following Jones and Launder [72], the diffusion term is modelled as 
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The last two terms on the RHS of Equation (41) are approximated as 
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Assuming the turbulent diffusion due to pressure fluctuations is small, 
Hanjalic [73] has found the following form may provide the basis for satis- 
factory predictions for high Reynolds number flows, 



(pUjE) 
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(44) 


The recommended empirical coefficients are 

= 0.09, c £l = 1.44, c £ 2 = 1.92, o t = 1.3, a k = 1.0 (45) 

This model is valid only in the turbulent flow regime where viscous effects 
are negligible. The wall function approach, that uses a form of the law of the 
wall to approximate the solution for the near wall region where viscous effects 
are not negligible can be used with the above model. This will be discussed in 
detail in section 2.2.3. 


2.2.2 The low-Reynolds number A-c model 

Jones and Launder [72] extended the two equation turbulence model to pre- 
dict the flow in the viscous sublayer. The model includes the viscous diffusion 
of A and e and the coefficients in Eqs. (38) and (44) are assumed to be de- 
pendent on turbulent Reynolds number. The approximate transport equations 
of Jones and Launder's low Reynolds number turbulence model(JL) [72] with 
some constants modified by Launder and Sharma [74] are summarized below. 
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The turbulent viscosity, n , , is given by 


kt= c ufaP — 


(48) 


where the dissipation rate, e, is equal to £ + 2v (dk ]l2 jdx k ) 2 , and the function 
f is introduced to produce the effect of molecular viscosity on the shear stress. 
The empirical coefficients adopted are 


fa = ex P 


-3.4 

(1 + Re i /50) 2 


f & = 1.0- 0.3 exp {-Re;) 

,2 


Re, = 
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(49) 


The above equations with the appropriate boundary conditions for h, and 
u, — k = £ = 0 at the wall comprise a complete model. 

Chieng and Launder [3] applied this low Reynolds number model to the 
flow in a sudden pipe expansion, and the predicted heat transfer rates in the 
vicinity of the reattachment point were too high by up to a factor of 5 as a 
result of too large a length scale near the wall. Yap [1 1] added a source term 
S £ to the right hand side of the transport equation for e 
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where y is the distance from the wall, and c t — 2.5. Yap found that inclusion 
of this source term gave Nusselt number predictions comparable with the ex- 
periment and improved the Reynolds number dependence. This term was also 
used in the present study. 


2 . 2.3 Wall functions 

For many turbulent wall boundary layers the inner portion of the flow has 
the logarithmic law of the wall behavior. With the wall function approach, the 
equations are solved on a relatively coarse grid, and the near wall region is re- 
placed with formulas based on the logarithmic velocity and temperature pro- 
file. The logarithmic velocity profile is given as 

u + = -L In (Ey) (51) 

where k is the von Karman's constant( = 0.4l), and E is equal to 9.0 for a 
smooth wall. 

The skin friction coefficient is defined as 


C L 

2 




(52) 


Substitution of Eq.(51) into the Equation (52) with the definition of local 
Reynolds number 
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Re = 


uy 
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(53) 


yields a relation for momentum transfer to the wall. 
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(54) 


The logarithmic temperature profile is given as 

7^ = Pr t {u + + P) (55) 

where P is an empirical function of the molecular Prandtl number. Substi- 
tution of the above equation into the definition of Stanton number 


St 




puc p {T w - T) 
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u + T + 


yields the relationship for the heat flux to the wall as 



(56) 


(57) 


For a smooth wall, Launder and Spalding [2] proposed a simpler version of 
Jayatilleke's correlation [75] of the P-function, 


P = 



(58) 


Convection and diffusion of turbulence energy are found to be nearly always 
negligible in the vicinity of a wall [76]. The balance between the production 
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and dissipation terms in the turbulence kinetic energy equation yields the 
boundary values of k and £ as 


k = 





£ = 


Jl± 
Q y 


(59) 


(60) 


where C d is a constant equal to 2.55. 

The above wall functions have been tested for channel flow in the present 
study and satisfactory results were obtained for both momentum and heat 
transfer. They were used for flow in a sudden pipe expansion to supply fully 
developed conditions at the inlet. 

Although the above wall functions correctly predict zero shear at rcattach- 
ment points, they also predict zero values of heat flux and turbulence kinetic 
energy. However, experimental measurements show that these quantities have 
maximum values at reattachment points. This deficiency can be removed by 
choosing kj, 12 as the velocity scale, which was first proposed by Launder and 
Spalding [2], This is discussed in the next section. 

Wall functions have been derived above for smooth walls. For rough walls, 
only the empirical functions need be replaced. For example, Han's [27] 
roughness functions for transverse rectangular ribs reduce to 


E = 


exp[0.39 (L / e)° 53 ] 


(61) 
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( 62 ) 


0.28 I n ci 

P = 5.05 e + - 0.97 


2.2.4 Improvements in wall functions 

As pointed out earlier, the wall functions based on the velocity scale of 
\/t Jp are not suitable for recirculating flows because of incorrect predictions 
of heat transfer rates and turbulent kinetic energy at reattachment points. The 
problem can be avoided by adopting the proposal of Launder and Spaldifig 
[2] in which the flow in the near wall region is assumed to be in local equilib- 
rium. The logarithmic velocity and temperature laws are used with the non- 
dimensional parameters replaced with 
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In the turbulent kenitic energy equation, the dissipation term used is the aver- 
age over the control volume near the wall. 
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For the equation of dissipation rate e , the value is evaluated under local 
equilibrium condition as 
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(65) 


= Jl 4 P 

* «y P 

The above wall treatment has been tested in the present study for a 
backward-facing step problem, and the heat transfer rates in the neighborhood 
of reattachment points were poorly underpredicted (as has been shown by 
many investigators). 

Ciofalo and Collins [23] extended the Johnson and Laundcr's proposal 
[4] of varying viscous sublayer thickness with the level of turbulent kinetic 
energy diffusing into the sublayer. Their approach is to relate the non- 
dimensional viscous sublayer thickness to the near-wall turbulence intensity 
rather than to the near-wall profile of the turbulent kinetic energy as in the 
proposal of Johnson and Launder. This proposal is implemented as follows. 
In the equilibrium boundary layer, the nondimensional turbulent kinetic en- 
ergy, k + = k I u} , can be expressed as 

k+ - c;' l2 (Kr) y + ^yl ( 6 «) 

y^o 

k* = c" 1 ' 2 y + >y+ (67) 

where y+ is the nondimensional viscous sublayer thickness taken as 1 1 .225. 

Using the above profiles of k + and those of u+ in the viscous sublayer and 
in the fully turbulent region, the turbulence intensity in the equilibrium 
boundary layer can be expressed as 
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( 68 ) 



With the calculated sublayer thicknesses, _y+ and j 4 , the constant E and the 
P function in the logarithmic laws can be reevaluated. 


35 



I 


2.2.5 Wall function method for variable properties flow 

Viegas and Rubesin [77] extended Chieng and Launder s wall function 
approach [3] to include the effect of compressibility for flow over adiabatic 
surface. Supersonic flow over a flat plate and two cases of normal shock-wave 
turbulent boundary layer interaction at transonic speed with and without sep- 
aration were computed with k — e model of Jones and Launder [72] and of 
Chien [78], and k - co 2 model of Wilcox and Rubesin [79 ]. Computations 
were performed with both the wall function method and integration to the 
surface with all three models, and comparison with experiments were generally 
better with the wall function method, provided that the first two mesh points 
lie between the buffer layer and the wake. 

Viegas et al. [80] improved the wall function approach to nonabiabatic 
conditions and relaxed the criteria for the placement of the near-wall mesh 
points so that they can lie in the viscous sublayer. Computations for attached 
and separated flows over both adiabatic and nonadiabatic surfaces with Mach 
numbers, 0.875 < Ma <2.85, gave good agreement with the experimental data. 
However, they point out that the nonadiabatic contribution to the wall func- 
tion had a small effect on the local temperature due to small surface heat fluxes 
and that the wall function can be in considerable error. They comment that 
improvements in shear and eddy viscosity modeling within the first mesh vol- 
ume will be required for consideration of very cold wall cases. 

Following Chieng and Launder's approach, Viegas et al. assumes that the 
turbulence Reynolds number, Re v = yjcl l2 jv , a universal constant equal to 20 
which would require a cubic equation for k v . This assumption has led to poor 
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prediction of heat transfer rates (20-30% lower than the experimental data) in 
abrupt pipe expansion flows [4J. In this study, the approach of Ciofalo and 
Collins [23] is used to calculate the viscous sublayer thickness, and the as- 
sumption of Re v = 20 is abandoned. However, the calculation of P function 
from the thermal sublayer thickness is not used, and the surface heat flux is 
calculated from a reduced energy equation as will be shown later. The effects 
of compressibility was added following the approach of Viegas and Rubcsin. 

The conventional law of the wall is extended to compressible flow with the 
use of the van Driest transformation, and the effective u + is 



and the effective kinetic energy of turbulence is shown to have the density 
scaling as [79] 


k rr= 

‘ Pw 

The nondimensional y coordinate is chosen as 

+ y * t 
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(74) 


(75) 


and the friction velocity is defined as 
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With the above relations, the mean velocity profile is given by 
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and a mean value of z was employed as 



T + T 
1 V “ L W 


(80) 
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In the fully turbulent region, neglecting the molecular diffusion of heat leads 
to the following temperature distribution. 
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(82) 


T = T w - ^ [ Pr t (-j-), (u - u v ) + Pr (4-)/ « v ] 


where 
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and a mean value of t is taken as 
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Substitution of Eq. (79) and Eq. (82) into Eq. (78) produces the following re- 
lationship for the effective velocity. 
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The velocity at the viscous sublayer is obtained from the momentum equation 
at the near-wall mesh as 
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If the near- wall grid point lies within the viscous sublayer, the wall shear stress 
is calculated using Eq. (79) and the assumption, /r//*w= T l T w 
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Given T w , the wall heat flux can be obtained from Eq. (82), and iteration is 
required among w v , T v , t w , and q w . 

Next, the mean generation and dissipation rates of k are derived as pro- 
posed by Chieng and Launder [3]. A simple expression for the mean dissi- 
pation rate I is obtained by neglecting the variation of k across the near-v;atl 
grid, but Chieng and Launder found that the prediction of local heat transfer 
rates were better with the variation of k taken into account. The turbulent 
shear stress is assumed zero in the viscous sublayer and undergoes an abrupt 
increase at the edge of the sublayer, varying linearly over the remainder of the 
near- wall grid. The mean generation of kinetic energy becomes 


P = 



( 88 ) 


w'here the subscript 1 denotes the north edge of the near- wall grid. Using the 
mean velocity distribution of Eq. (77) the integration leads to 
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For the evaluation of the mean dissipation, a linear variation of k is assumed 
between the edge of viscous sublayer and the second grid next to the wall, and 
a parabolic variation of k is assumed within the viscous sublayer as follows. 
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The dissipation rate of k is not zero in the viscous sublayer unlike the gener- 
ation. Use of the parablic profile of Eq. (90) gives 


e = 


2v vv k v 
Tv 


(92) 


where v in the sublayer is assumed to be a constant and equal to its wall value. 
In the fully turbulent region, e is evaluated as 


£ -- 


c iy 


(93) 


where c t — k/c^ 4 = 2.55 

The mean dissipation rate of turbulent kinetic energy in the near-wall grid 
is evaluated as 



(94) 


Upon evaluating the integrals, the following expression is obtained. 


+ _!_(iUj>/2 _ ft3/2) 2 a(k ll2_ k U2 ) +A (95) 
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where 
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a < 0 


Dissipation at the near-wall node is approximated under local- equilibrium 
condition as 


.3/2 




(97) 


2 . 2.6 Algebraic stress model 

The algebraic form of Reynolds stress equation can be obtained by neg- 
lecting gradients of 11$ Ik, which is true when jk is constant, and approxi- 
mately true when « ( w, jk varies slowly across the flow field. The resultant 
equation proposed by Rodi [81] assumes that the convective and diffusive 
transport of ujUj are analogous to those of turbulent kinetic energy as follows. 


Ui 


dujUj 

dxb 


D U = 


u i u j j, dk 
k k dxi 


U t Uj 


D(k) 


(98) 


where Dy is the diffusive transport tensor and D(k) is the diffusion of k. The 
transport equation for k and e are modelled following Daly and Harlow [82], 
as follows: 
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Combining the equations (98) and (99) with the Reynolds stress equation, the 
algebraic equation for UjUj is obtained as 


UjUj 

T” 


(P-e) 


P ij-'T 6 ij £ + ^ij 


(1U1) 


where P is the production of k, and 4> i} is called the pressure-strain tensor. The 
generation tensor of ufy is given by 



- u i u k 



- UjU k 


dUj 

dx k 


( 102 ) 


The pressure strain term is modelled into two component, <j> Uti which involves 
only turbulence quantities, and which involves products of turbulence 
quantities and mean rates of strain. The component has been refered to 
as the "return-to-isotropy" term, and is modelled by Rotta [83] as 

(l03) 


The mean strain part of the pressure strain term, (f> ij2 has been modelled by 
Launder et al. [84] as 


-'2 (/>,;-} hjfj 


(104) 
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where c, and c 2 are 1.8 and 0.6 respectively. Shir [85] proposed that the term 
<f> y] be modified to account for the near-wall effects by adding to it the fol- 
lowing correction: 

<f>ij, i = ci n k n m b ij~\ Wi n k n j - y Wj n k ( 1 ' 05 > 

where n t is the unit vector normal to the wall, r,- is the distance from the wall, 
and / is the length scale, k 3l2 lc,e. The /function acts as to diminish the influ- 
ence of the near-wall correction as the distance r increases. Shir's idea v.’as 
extended by Gibson and Launder [ 86 ] to model a near-wall correction for the 
term <y >2 . 

<hj ,2 = c 2 n k n m S ij ~ \ <t> ik,2 n k n j ~ \ <t>jk,2 n k "/)/( n/ 77 ) ( 106) 


The final algebraic equation for the Reynolds stresses is 


u i u j — 3 v tj 


0 " C2)(t 


P+{C X - 1 ) £ 


(107) 


A proposal was made by Launder [87] to produce better prediction of turbu- 
lent shear stress in a free jet. His proposed equation is 


«/ u j 


T 5i J k+ jr 


(1 - c 2 ) 


p„-jL§ .. p] + — I — 

* 1 ” ) {l- C 2 ) 


1)(1 + a) + (a - /?) A + cj £ 

(&!// + tifj) 


(108) 


where 
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X = Pje and A = D{k)j e 


(109) 


a and /? arc empirical constants with values of 0.3 and -0.8 respectively. 


2.2.7 Closure 

The wall function approaches were tested in the present study for abrupt 
pipe expansion flow and flow downstream of a backward-facing step, and 
proved to be inadequate(especially when the step height is too small to locate 
sufficient number of grid points). The low-Reynolds number model with 
modification by Yap has received most attention in the present study. 

It has been known that the Boussinesq-viscosity hypothesis cannot simulate 
the level of anisotropy of normal stresses resulting from curvature in flow. A 
related example is the turbulence-driven secondary motions causing bulging 
of the velocity contours towards the corners in straight, non-circular ducts and 
open channel. The motions has been known to have a pronounced effect on 
the shear stress and heat transfer in the corner region and cannot be predicted 
by an isotropic eddy-viscosity model. Numerical studies are currently under 
progress w r ith the algebraic stress model. 
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Chapter III 

NUMERICAL SOLUTION 


3.1 PHOENICS 

The present study uses PHOENICS(Parabolic Hyperbolic Or Elliptic Nu- 
merical Integration Code Scries) developed by CHAM Ltd. It is a general 
purpose computer program for the analysis of fluid-flow, heat transfer, 
chemical-reaction and related phenomena. PHOENICS can solve single-phase 
or two-phase parabolic and elliptic problems in cartesian, cylindrical polar, 
and curvilinear body-fitted coordinate systems. Dependent variables in the 
governing equations are allowed to vary in one, two or three dimensions and 
in time. 

PHOENICS consists of two main computer codes and an auxiliary code. 
The main codes are a pre-processor called SATELLITE and a processor called 
EARTH. SATELLITE is an interpreter that converts instructions provided 
by the user into a data file for EARTH. EARTH is the main flow-simulating 
software that executes the corresponding computations. Various data-setting 
can be made by the user in GROUND which is a subroutine of EARTH. The 
auxiliary code is called PHOTON. It is a graphics program that presents the 
computed grid and flow pattern on the screen. 
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3.2 SOLUTION ALGORITHM 


3.2.1 Grid configuration 

In PHOENICS, staggered grids are used. Pressure and other scalar de- 
pendent variables including the turbulence kinetic energy and dissipation rate 
are stored at the cell nodes whereas velocities are located at the cell faces. For 
U velocities, the control volumes are displaced in the positive x direction by 
one half cell so that, for U p , the west boundary of the cell passes through the 
node point P, and the east boundary passes through the node point E. Grid 
node locations and the staggered grids for velocities are shown in Figures 3.1 
and 3.2. The V velocity control volume is similarly displaced by one half cell 
in the y direction. 

The grid nodes are placed in the center of the control volume, rather than 
the control volume faces being placed half way between grid points. This grid 
arrangement is convenient with the solution domain containing porosities, 
since the control volume faces can be located along the boundaries of the 
blockage, whereas the latter grid arrangement would require setting up the grid 
nodes first and cell boundaries to coinside with the boundaries of blockage. 

The governing equations are integrated over individual control volume or 
grid cell to arrive at the discretised equations. 


47 



3.2.2 Discretisation of the governing equations 

PHOENICS provides solutions to the discretized versions of sets of differ- 
ential equations having the general form [88] : 

-fkrjPiM+V ■(nPi V<h-nV 4 ,F<h)=r,S, ( 110 ) 

at 

where, 

t stands for time 

r i stands for volume fraction of phase i 
stands for density of phase i 

(f>i stands for any conserved property of phase i , such as enthalpy, 
velocity, mass fraction of a chemical species, etc. 

V stands for velocity vector 

T H stands for the exchange coefficient of in phase i 
Si stands for the source rate of <j>i 

Integration of Equation (1 10) for a single phase over the whole volume of the 
domain of surface area A, followed by application of the divergence theorem 
yields: 

J = jjjv v (lll) 

Rewriting the convective and diffusive fluxes normal to the surface area at 
e, w, n, s (east, west, north, and south of the control volume faces), the 
equation in two-dimensions becomes: 
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J e K ^ ^ n J s 


V v 


( 112 ) 


with 
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J e = (pU) e (^ D fiy ns) 4> e ~ ^ <f> ,e p ns) q x U 


vv 


e \'p ^ /tr' t - e 

J w = (puu&yns)* 
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i d(j) 

r 4>,w( r p 5 yns)^\w 
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r 4>,nri 5x ew)-fy\n 

i d<t> 

^(ppc ^ x ew) b 


( 113 ) 


and Sy ns and Sx ew indicate the distance between cell volume walls, n and 5 , and 
the distance between cell volume walls, e and vv respectively, as shown in Fig- 
ure 3.1. j = 0 and 1 corresponds to cartesian and cylindrical coordinates re- 
spectively. For the convective and diffusive flux at the west face, J e , the 
transport properties at e can be evaluated cither by arithmetic averages of 
those on either side of the cell faces, or by harmonic mean averaging. In the 
present study, arithmetic averages have been used to approximate the trans- 
port properties. 

The value of <f> is assumed to vary linearly between grid nodes in the ap- 
proximation of d(pldx \ e , and the use of central difference scheme results in: 


dcf> 

dx 


(</>£— 4>p) 


5x 


( 114 ) 


EP 


In the evaluation of (f> e in the convective flux term in J e , there are many alter- 
natives. In the upwind scheme, a stepwise variation of 4> is used, and the flux 
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of variable <f> across the cast cell face is taken as the product of the mass flux 
and the value of <fr at the upwind node. Thus the flux is: 


p p U e Sy n s4> P for i/ e >0 

PE U eh ; ns ( i ) E for U e<° 

In the hybrid scheme that was developed by Spalding [89] , the value of (j> e is 
evaluated by the central difference scheme when the cell Peclet number, Pe, is 
less than 2 as follows. 


<l>e 


4 > p + 4 >e 
2 


(p U) p Sx pc 
I p ' — I < 2 
1 <t>,e 


(116) 


de Vahl Davis and Mallinson [90] have shown that the false diffusion coeffi- 
cient, arising from the use of upwind differencing scheme when flow cuts 
across grid lines at an angle, is given approximately by 

pUAx Ay sin 26 

ry = — — r - < ll7) 

4 ( A_y sin p + Ax cos 6) 


where 6 is the angle made by the velocity vector with the x direction. The false 
diffusion is the largest where the velocity vector is at an angle of 45° with the 
grid lines. 

A better approximation of the exact solution using a quadratic interpolation 
equations called QUICK was proposed by Leonard [91]. For a uniform grid 
spacing, the resulting formula is 

4> e = — - — - ^P + ^e) f° r 

2 8 (118) 

4 >e= - 2 — ■£" ^P~ ^£ + ^Ee) for t/ e >0 
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where the second term is like a correction to a linear interpolation proportional 
to the upstream-weighted curvature. EE is the cast neighbor grid point of the 
grid E. PHOENICS only has the hybrid scheme and the upwind scheme as 
options, and the present results were obtained with the hybrid scheme. The 
algebraic stress model with the QUICK scheme is currently under investi- 
gation. 

Substitution of the approximations for the profile of (f> at the cell face into 
the Equation (112) yields the set of equation of the following form. 

a p4>p = +jjv v (ll9) 

where i = E , IV, S, N for the hybrid differencing scheme. For scalar variables, 
the formula for the coefficients become 

a E = max (0, e — - a | p e U e by m | ) + max (0, - p e U e by m ) 

dx PE 

p 5 

a.yy = max (0, — a | p w U w by ns | ) + max (0, p w U w by ns ) 

dx WP (12 0) 

a N = max (0, | Xgw - a | p n V n bx^ | ) + max (0, - p n V n bx w ) 

fypN 

r 

ao = max (0, — - — — - a \ p s V s bx ew \) + max {0,p s V s bx^) 

fysp 

For a = 0, diffusive effects contribute irrespective of the value of cell Peclct 
number, and the upwind difference scheme is obtained. The hybrid difference 
scheme corresponds to a = 0.5. In the variable properties flow, the cell-face 
densities are evaluated using the upwing convention. Thus 
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P e = Pp for U e > 0 

Pe= PE for U e<° 

The expression used for the coefficient a E of U - velocity equation is 
ajr = max (0, d E - a | m E | ) + max (0, — m E ) 


(121) 


( 122 ) 


where 


m E = 


Pe ns ( ^ e + ^eF) 


d E ~ 


^ E ns 
& x eeE 


(123) 


eE is the east cell face of the grid point E, and other coefficients have the 
similar form. 

In the present study, the sets of linear equations are solved by the 
TDMA(tridiagonal matrix algorithm) which is described in section 3.5. 


3.2.3 Source term linearization and boundary condition 

The source term needs to be a linear function of </> in order to have the 
whole discretised equation in the linear form. The nodal value of source terms 
are supposed to prevail over the whole of the cell volume, and the source term 
over a control volume can be written as 

Jjvv = S c +C,^ (124) 
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The proper linearization is to have the coefficient C ^ be less than or equal to 
zero in order to assure that the coefficient a P remain positive to avoid diver- 
gence of solution iterations. For the momentum equations, the pressure gra- 
dient is added to the source term, and all the boundary conditions also enter 
the discretised equations by way of the source terms. 

The boundary condition on the wall for the momentum equations enters the 
source term as the product of wall shear stress and the surface area of a wall 
bounded cell face. For the laminar flow problem or a lovv-Reynolds number 
turbulence models, the wall shear stress is simply the product of the near wall 
velocity gradient and the fluid viscosity. In the case of wall functions ap- 
proach, the wall shear stress can be evaluated from the law of the wall as in 
section 2.2.3. 

The constant wall temperature boundary condition enters the source term 
similarly as the product of temperature gradient and the fluid thermal 
conductivity near the wall. For the constant wall heat flux boundary condi- 
tion, the coefficient Q, is set to a small number, 10 -10 , to prescribe a fixed-flux 
source term. 

In PHOENICS, the source term is always expressed as a linear function of 
the dependent variable 4> as [92] 

S+ = ~ Mt (125) 

where C ^ is a link coefficient which relates the source term to the difference 
between in-cell value of 4> and the boundary value, and V ^ is the boundary 
value for variable (p. f t represents some geometric factor such as area of the 
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computation cell. This equation is appropriate for a boundary across which 
there is no mass flow, and source or sink terms only arise due to diffusion. 

For an inflow boundary, the mass source is expressed as a linear function 
of pressure: 

m=C p {V p -p)f t (126) 

The source of any other dependent variable 4> is 

Sff, = thV^ + ~ 4>)ft ( 1 ^) 

where the first term represents the amount of </> convected into the computa- 
tion domain, and the second term represents the source of <j> due to diffusion. 

= 0 is appropriate if diffusion effects are negligible. In order to prescribe 
a flux boundary condition with Eq. (125) C ^ is set to a small number, 10 10 , 
and Vq is multiplied by 10 +, ° . 


3.2.4 Solution of hydrodynamic equations 

There exists various methods of treating the pressure-velocity coupling be- 
tween the mass and' momentum conservation equations. PHOENICS uses a 
variant of the SIMPLE algorithm [93] called SIMPLEST. The major differ- 
ence between the SIMPLEST and SIMPLE algorithm is that in the former the 
coefficients for the momentum equations contain only diffusion contributions, 
and the convection terms are added to the linerized source terms. This implies 
that, in the absence of diffusion, the momentum equations are solved by a 
Jacobi point-by-point procedure instead of the line-by-line procedure [94]. 
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The essential idea is to use the continuity equation to derive an equation for 
pressure correction to be added to the current iteration value of pressure which 
will tend toward the correct pressure and also satisfy continuity equation. This 
is summarized below. 

The discretized momentum equation with the staggered grid for the velocity 
components can be written as 

a e u e = ^ a { iq + b + (p P — pp) A e ( 1 28) 

An incorrect pressure field in the momentum equation yields a velocity field 
that will not satisfy the continuity equation. In order to correct the guessed 
pressure, the following steps of the SIMPLE algorithm are used, as proposed 
by Patankar and Spalding [93]. 

The correct pressure and velocity are assumed to be in the form: 

P = P* + P' 

u = u* + u' (129) 

v = v* + v' 

where the starred variables are imperfect values and the primed variables are 
correction terms. Subtracting the discretized momentum equation based on 
incorrect pressure and velocity fields from Eq. (128) results in : 

a e u e = ^ a i u i + ~ 1 

The first term on the right hand side of Eq. (130) is neglected since the con- 
verged solution given by this algorithm does not contain any error resulting 
from its omission. The resulting velocity correction formula is 
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( 131 ) 


< = ^(PP-PE) 

e 

Equation (131) can also be written as 

u e = u *e + ~lT~(Pp - Pe) ( 132 ) 

e 

Equation (132) can be used to correct velocity from the pressure corrections. 
Substitution of Eq. (132) and the corresponding equation for the velocity 
component, v„, into the discretized continuity equation of the following form . 

[ ( pu ) e - ( pu ) w ] by - [ ( pv )„ - ( pv ) s ] 5x = 0 (133) 

yields the pressure correction equation. 

a p Pp = ^ a iP{ + & (134) 

where 

b' = l(pu\- {pu*) e ~] by + t(pv*) s - (pv*)J 5x (135) 


b’ is a mass imbalance in the control volume due to the fact that the current 
velocities do not satisfy continuity. 

The overall solution procedure consists of the following steps. 

1. Guess the pressure field. 

2. Solve the momentum equations to obtain the velocity field based on the 
guessed pressure field. 

3. Solve the pressure correction equation. 
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4. Solve for the new velocity field using the velocity correction formula. 

5. Solve the discretized equation for other variables ( temperature and 
turbulence quantities ). 

6. With the new pressure field, return to step 2 and repeat the procedure 
until a converged solution is obtained. 

3 . 2.5 Tridiagonal matrix algorithm 

The tridiagonal matrix algorithm can be used to solve any set of equations 
when the matrix of the coefficients of the equations consists of nonzero coeffi- 
cients aligned along three diagonals of the matrix. The discretised equation 
(119) can be rearranged in the following form. 

a j = b j fy - 1 + c ] ^j + 1 + d j ( l36 ) 

where the subscript j refers to the grid node P and varies along a chosen line. 
The recurrence relations are 


C, = 


D; = 


9 

a j~ b j C j- 1 
dj+bjDj-X 

a r b j c j - 1 


(137) 


and the second stage consists of a back-substitution by 


$j ~ c j ^j + 1 Dj 


(138) 


For the first equation 
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(139) 


and 

/ 

<t> N = D n (140) 

A sweep along the line of nodes in the transverse direction is performed, and 
the calculation is performed on the next parallel line of nodes with updated 
values of 0. Upon completion of sweeps in the transverse direction, a simUar 
procedure may be performed in the sweep direction. 


3.3 COMPUTATIONAL ASPECTS 
3.3.1 Convergence criterion 

The criterion of convengence of the numerical solution is based on the ab- 
solute normalized residuals of the equations that was summed for all cells in 
the computation domain. The mass residual, or the imbalance of the conti- 
nuity equation, is defined as 


RES 



(141) 


where b' is defined in the Equation (135). The residuals associated with other 
dependent variables are defined as 
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RES = 


(142) 


S'Z a i 0/ + S ~ a p4 > p\ 

where, for the momentum equations, M in ^ is the total inflow of momentum 
and for the equations of the turbulence kinetic energy and dissipation rate, it 
is the product of total volumetric inflow and the inlet value of k and e respec- 
tively. The solutions are regarded as converged when these normalized resi- 
duals become less than 10~ 3 for the continuity equation and 10~ 2 for other 
variables. 

In the case of the turbulent flow in a channel with ribs, a typical output had 
the normalized absolute residuals of le — 5, le—2 , l.5e— 3, 1.67, 1.67, 
4.4e— 3 , for continuity, V momentum, U momentum, k, e, and energy 
equation respectively, after 1500 sweeps. 

In addition to the whole-field residuals, the average friction factor and the 
average heat transfer coefficient were monitored at every 50 sweeps for the 
problem of flow in channel with ribs. The relative errors of the average friction 
factor and the average heat transfer coefficient defined below was less than 1% 
for all cases studied. 

j n _ 

rel. error = ™ % (143) 

4> n 

In the case of flow in an abrupt pipe expansion, flow over a backward- 
facing step, and the impingement cooling, the local heat transfer coefficients 
were also monitored at every 100 sweeps. When the solutions were well con- 
verged, the computed local heat transfer coefficients were almost invariant. 
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For example, the change in the stagnation heat transfer coefficient for the 
impingement cooling problem was less than 1% over 1000 sweeps once the 
heat transfer coefficient profile was converged. 

For the flow in a tube with transverse rectangular ribs, an increase in the 
number of grids in the axial direction of 30% changed /and St by about 4% 
and 3% respectively, and improved agreement with experiment. For the most 
of the computation domain, at least 3 grid nodes were placed within the 
viscous sublayer, and the grid was stretched in x or y direction by a factor of 
1.1 -1.5 carefully over various parts of solution domain. The effect of further 
grid refinement near the wall was negligible for the prediction of flow field. 
For the computations with wall function methods, most of the near-wall grids 
were placed in the fully turbulent region, y + > 30, and extensive grid refine- 
ment tests were not performed. 


3.3.2 Under-relaxation 

The present computation involves equations that are nonlinear and strongly 
coupled, and under-relaxation is required to achieve the overall convergence 
to a solution. In the present calculations under-relaxation was applied in the 
following manner involving the false time step 5t f . 

{ ap+ '^) <1>p = Y a ^> +s+ i^ <t>p (144) 

where (p* P is the previous sweep value of (f> at point P. If <5/ is small the terms 
containing it will be large and tend to dominate the equation, implying 
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<3 i>p = 4>*p 


(145) 


Thus, the smaller the value of bt f the less <p can change from sweep to sweep. 

An alternative form of under-relaxation called linear under- relaxation was 
used for pressure as 

p = p x + a {p — p x ) (146) 

The turbulent viscosity was also under-relaxed in the similar manner. 

p t = p* t + a{p t -p*) (147) 


3.3.3 Computational time 

For the flow in the abrupt pipe expansion of djD = 0.8 and Re = 20,130, 
a grid of 71x37 for the low- Reynolds number turbulence model required about 
620 seconds on the IBM ES/9000 model 900 in order to achieve the conver- 
gence criteria. The RES of the turbulence kinetic energy and the dissipation 
rate were not lowered below 0.5, but the local Nusselt number prediction was 
almost invariant, and the iteration was terminated at 1500 sweeps. A grid of 
59x21 for the wall function methods and the same problem required less than 
half the time required for the low-Reynolds number model to achieve approx- 
imately the same convergence criteria. 

In the case of channel with ribs, the average friction factor and the average 
Stanton number were monitored, and a typical run with a grid of 96x66 re- 
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quired about 1300 seconds. About the same amount of computing time was 
required for the jet impingement cooling problem with a grid of 58x74. 


3.4 COMMENTS ON PHOENICS CODES 

PHOENICS has many capabilities for simulation of various fluid-flow, heat 
transfer, and related phenomena. For example, curvilinear body-fitted coor- 
dinates option can be used to generate grids for complex geometry. However, 
the code is weak in turbulence modelling; the standard k — e model is the most 
advanced one provided. For the low Reynolds number k - e model, the addi- 
tional terms in the transport equations of k and e have to be coded as a source 
term. The two-dimensional algebraic stress model formulation requires as 
many as ten extra source terms to be coded in each of the momentum 
equations. Due to inaccessibility of the source code and its large size, modifi- 
cation can become very difficult. The advantage of PHOENICS over other 
specific codes may be its capability to solve diverse problems, but due to its 
large applications the code is difficult to debug. In addition, the upwind 
scheme and hybrid scheme are the only differencing schemes available in the 
current version of PHOENICS(1.4). Thus, the code may also suffer from lack 
of solution accuracy due to numerical diffusion in many problems. 
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3.5 BOUNDARY CONDITIONS 


3.5.1 Flow in a sudden pipe expansion 

Inlet conditions were fully developed profiles for k, £, and the axial velocity, 
U, obtained from the straight pipe flow solutions of the standard k-z model 
with wall function boundary conditions of Eqs. (54), (57), (59), and (60). The 
inlet condition for enthalpy was a uniform profile. The outlet was located 
about 200/?( = 60£>) downstream from the step so that its influence on the main 
flow would be negligibly small. The boundary conditions imposed were 
outlet : 

P = 0 

= 0 for f = u, v, k , £ and T 
dx 

symmetry : 

- = 0 for (f> = u,v,k, z and T 
dr 

wall : 

u, v, k, £ = 0 
constant wall heat flux 

The wall boundary conditions were an adiabatic wall on the entry pipe 
( - 0.5 D < x < 0) and on the side wall of the step, and constant wall heat flux 
on the expansion tube as in Baughn's [9] experimental setup. 
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3.5.2 Flow over a backward-facing step 


The same boundary conditions for the flow in a sudden pipe expansion ap- 
ply to this problem with cylindrical coordinates replaced by cartisian coordi- 
nates. The inlet profiles for k , £, and the axial velocity U were obtained from 
the two-dimension channel flow solutions of the standard k — e model with 
wall function boundary conditions. 


3.5.3 Flow in channel with rectangular ribs 

The fully developed condition of du/dx = 0 and v = 0 cannot be applied to this 
problem since u varies continuously with x and v is not zero. But, the flow 
field will repeat itself in a succession of cross sections that are separated by the 
pitch length, L, as shown in Figure 3.3, sufficiently far downstream. One 
possible approach is to use periodicity boundary conditions as proposed by 
Patankar et al. [31]. Their approach is summarized below. 

The velocity components are assumed to behave periodically as 


u,{x,y) = ufx + L,y) = u t {x+2L,y) = (148) 

The periodicity condition for pressure is 

p{x,y ) — p(x + L,y) = p{x + L,y) — p(x + 2 L,y) = (149) 

and the pressure is assumed to be composed of 

p(x,y) = - x + P(x,y) (150) 

where /? is defined as 
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(151) 


_ P(x,y)-p(x+ L,y) 

P ~ L 

The fix term is related to the global mass flow, and P{x,y) is related to the 


detailed local motion, and is assumed to be periodic as 

P{x,y) = P{x+L,y)= P{x + 2L,y) = (152) 

The fully developed temperature profile condition of dTjdx = 0 for con- 


stant wall heat flux boundary condition cannot be applied to the present 
problem for two reasons. The first is the non-uniform heat transfer surface 
area which precludes uniform heat addition to the fluid, and the second is the 
nonzero axial conduction term, d 2 Tldx 2 . The temperature field is assumed to 
be composed of a component due to the heat flux plus a periodic component. 

T(x,y ) = yx + T{x,y) (153) 


where y is defined as 

7fx+ L,y)- T\x,y) = Q 
^ ~~ L mc p L 

A 

and T is periodic as 

T(x, y) = T( x + L,y)= T(x + 2 L,y) = . 


(154) 


(155) 


Substitution of the Eq. (150) and Eq. (153) into the momentum and energy 
equation yields a fi term on the right hand side of the x-momentum equation 
and a - uy term on the right hand side of the energy equation respectively. 
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The pressure gradient /? generates a corresponding mass flow, and its first 


value is assigned as 


P 


2 

P u b 

2D 


(156) 


where /is obtained from the flow in a pipe with repeated- rib roughness [26]. 
The value of /? can be iterated upon until the solution converges, as proposed 
by Lee [30]. 



1 +<x 




(157) 


where the subscript o refers to the value of the previous iteration, and a is an 
overrelaxation factor. 

Patankar et al. have used this approach for the laminar flow problem, but 
with at least two more equations(turbulent kinetic energy and dissipation rate) 
in the present problem, further complication of the convergence problem was 
avoided, and a simpler approach to the periodic boundary condition was used. 
A fixed number of inner iterations were performed for given inlet conditions, 
and the calculated outlet values of velocities, enthalpy, kinetic energy, and 
dissipation rates were substituted as inlet conditions for the next outer iter- 
ation. A 1 /7th power law profile was given for the axial velocity for the entire 
field as an initial guess in order to accelerate convergence to a fully developed 
condition. For this approach, the computation domain in Figure 3.3 was 
modified such that the inlet was located at six slabs before the right end of the 
first rib for the flow in a channel. The values at slab NX-1 were substituted 
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as the inlet conditions after each outer iteration. For the flow between parallel 
plates with wider ribs, the inlet was located 6 slabs before the right end of the 
first rib. The following boundary conditions were imposed, 
outlet : 

P = o 

- 0 for (f> = u,v,k, e and T 
dx 

symmetry : 

- = 0 for f = u,v,k, s and T 

dy 

wall : 

u, v, k, e = 0 

constant wall heat flux or constant wall temperature 


3.5.4 Jet impingement cooling 

For the present numerical computations, the following boundary conditions 
were used, 
outlet : 

p = 0 

= 0 for (j) = u,v,k, tandT 
dx 

symmetry : 

= 0 for f = w, v, k, e and T 
dr 

wall : 

w, v, /c, £ = 0 
constant wall heat flux 
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At the jet inlet, uniform profiles of turbulent energy, and energy dissipation 
rates were 

k in = iul (158) 

E in = k^KXdj) (159) 

where i is the turbulence intensity and X is the length scale constant. Typical 
values used were 0.5% and 30% for i and X respectively. Fully developed 
profiles of velocity, k, and £ obtained from pipe solutions were also used for the 
inlet condition. A schematic diagram of the computation domain is shown in 
Figure 3.4. 


3.6 COMPUTATION OF FRICTION FACTOR AND STANTON 
NUMBER 

For the flow in a channel with rectangular ribs, The average friction factor 
is calculated from the pressure drop over one pitch length, 


/= 


P u l L i D h 


The local Stanton number is defined as 


St = 


lW 


P U b C p( T w ~ T b ) 


(160) 


(161) 


where 



(162) 


n= 



T udy 



The bulk temperature calculated by Eq. (162) was compared to that obtained 
from an energy balance in the abrupt pipe expansion flow, and they were al- 
most identical. Along the front and rear faces of the ribs, T b is taken as an 
average of the values at the upstream and downstream slabs. 

An average Stanton number is calculated as 


St„ 


P u b c p( T w ~ T b ) 


(163) 


where the average value of (T w — T b ) was obtained as 

\(T W -T b )dx 

T w - T b = -5 (164) 

where L is the rib pitch. Use of this definition of St m for a uniformly heated 
surface may give good agreement with average Stanton numbers obtained with 
a uniform temperature boundary condition, as described by Mills [95 ]. The 
average Stanton number including the front and back faces of ribs was also 
calculated using the entire length of the heated surface for L in Eq. (164). The 
difference between the two average St was found to be less than 5% in the 
present computations. 
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3.7 VARIABLE PROPERTIES FLOW 

The assumption of constant fluid properties is not adequate for large heat 
fluxes into the fluid, since all the physical properties depend on temperature 
and pressure. Fluid properties for the numerical calculation can be entered 
as power law approximations. The properties for low-pressure air can be ap- 
proximated within 4% in the temperature range of about 273K < T < 1500K. 
approximated as follows [96] 



0.67 


= Hr~) 


T; 


in 


k in 


1 in 


0.805 



(-JT-) 
1 in 


0.095 




165) 


Due to decrease in density with temperature, there is continuous acceleration 
of the flow, and also another effect is the increase of the fluid resistance at the 
wall with heating that causes thickening of wall boundary layer, thus reducing 
the core flow area. These effects must be taken into account in evaluating the 
wall shear stress. Assuming static pressure is uniform across the flow section 
and treating the momentum flux as one-dimensionl, the wall shear stress be- 
comes 


D djp + pulll) 
4 dx 


(166) 


In the present computation, the friction induced by the momentum change was 
taken into account in calculating the average friction factor over one pitch 
length. 

The specific heat was taken as a constant, and average viscosity and tem- 
perature ratio over the length of a pitch was taken as the average of the inlet 


70 


and outlet values. A fixed number of inner iterations were performed with a 
constant wall heat flux boundary condition, and the calculated values of ve- 
locities and turbulence quantities at the outlet of the computation domain were 
substituted as inlet conditions for the next outer iteration. 

Supercritical hydrogen properties are strongly dependent on temperature 
and pressure. Convenient curve-fits have been supplied by Baek [97]. These 
curve-fits have been incorporated in PHOEN1CS and can be used in the future 
study of complex turbulent flow involving supercritical hydrogen. 
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Chapter IV 

RESULTS AND DISCUSSION 

4.1 FLOW IN A SUDDEN PIPE EXPANSION 

For the present computations, expansion ratios of 0.4 and 0.8 and a down- 
stream Reynolds number range from 10570 to 76080 were chosen to study the 
effect of Reynolds number and the expansion ratio on the performance of the 
turbulence models. The inlet conditions were fully developed profiles for all 
the variables obtained using wall function boundary conditions, except 
enthalpy, which was uniform. Air was chosen as the fluid, and properties were 
evaluated at 293 K. The wall heat flux used was 100W I m 2 which was in the 
range of Baughn's experiment [9]. The turbulent Prandtl number was fixed 
at 0.9. A typical grid used for djD — 0.4 was NY = 59 and NX = 71 with 75% 
of the grids located between the wall the the top of the step in the radial di- 
rection. 

Figure 4.1.1 shows the velocity profiles for the expansion ratio of 0.4 and 
0.8. The local Nusselt number distribution normalized by the Dittus- Boelter 
relation for djD = 0.4 and Re= 12310 computed by three methods, the low- Re 
model with modification by Yap, wall function method of Collins, and the 
standard wall function method are shown in Figure 4.1.2a The low-Re model 
predicts the overall Nusselt number distribution the best except for an under- 
prediction of reattachment length. The standard wall function and the method 
of Collins underpredict the maximum Nu at the reattachment point by about 
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40% and 23% respectively. The heat transfer rates in the recirculation zone 
are predicted better by the wall function methods, as have been reported by 
other investigators. 

Figure 4.1.2b shows that the corner eddy near the step can be predicted by 
the low- Re model, but the maximum wall shear stress in the recirculation re- 
gion is about 2 times larger than that by the methods of Collins. Experimental 
data for the wall shear was not available, and a comparision could not be 
made. The skin friction and the Nusselt number are both low near the st<;p. 
However the skin friction approaches a maximum about 5 step heights up- 
stream of reattachmcnt and declines to zero at reattachment while the Nusselt 
number reaches to a maximum near the reattachment. The skin friction and 
the heat transfer coefficient behave differently in a recirculating flow, and the 
Reynolds analogy does not hold in a reattaching or recirculating flow. Thus, 
the wall function approach of relating the flow with the wall shear stress fails 
in the neigborhood of reattachment. 

Figure 4.1.3 shows the predictions of local Nu distribution for Re = 23,210 
and 40,750. The low- Re model overpredicts the maximum Nu of the exper- 
iment data by about 25% for Re = 40,750, whereas the Collins method gives 
a good prediction at the reattachmcnt and in the recirculation region. How- 
ever, the predictions by wall function methods are consistently poor in the re- 
development region. The values of Nu db for Re= 12,310, 23,210, and 40,750 
are 37.6, 62.4, and 97.8 respectively. Figure 4.1.4 shows the computation re- 
sults by the low- Re model for Re = 10750 and djD = 0.8. For this expansion 
ratio, the step height is only 20% of the pipe radius whereas for djD = 0.4 the 
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step height is 60% of the pipe radius. Also, use of the wall function boundary 
conditions was not practical for this low Reynolds number, because not 
enough grid points could be placed below the step if the first grid point adja- 
cent to the wall was to be selected to lie within the fully turbulent part of 
boundary layer. The prediction of local Nusselt number shows a good 
aggreement with the experimental data over the entire field. Figure 4.1.5 
shows predictions by the three methods for Re = 20,130 and 39,300. The wall 
function methods perform poorly for this small step height and the Reynolds 
number range. Prediction for Re = 76,080 in Figure 4.1.6 shows that the 
low- Re model overpredicts the maximum Nu by about 18% and again indi- 
cates the need for correction at high Reynolds number. 


4.2 FLOW OVER A BACKWARD-FACING STEP 

Vogel and Eaton [19] have performed a detailed study of fluid flow and 
heat transfer for flow over a backward-facing step. Numerical computations 
have been performed to investigate the performance of the low- Re model and 
the wall-function methods. The experiment of Vogel and Eaton was performed 
had a development section of length 2.5m, and transpiration was used to vary 
the boundary layer thickness. For the present computations, inlet profiles for 
the computations were obtained from the solution of the standard k — e model 
at the end of a channel of the length of 2.5 m (16.5D) , and a boundary-layer 
thickness of \.\e was obtained, where e is the step height. The outlet was lo- 
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cated at 54e downstream from the step. The expansion ratio was 1.25 and the 
Reynolds number based on the step height was 28,000. 

Figure 4.2.1 shows the vector plot, and Figure 4.2.2 shows the mean veloc- 
ity profiles at various nondimensional streamwise coordinates, x * = 

x — x r I x r where x T is the reattachment length. Predictions by the low-Re 
model and the wall function method of Collins are almost identical and are 
generally in good agreement with the experimental data. Figure 4.2.3 shows 
the comparison with the experimental data for the mean velocity profiles in the 
near-wall region that were measured by the traversing pulsed wall probe with 
two orientations. The experimental data at x* = 0.4 indicate an almost vertical 
profile whereas the profile at x* = 0.33 in Figure 4.2.2 shows a gradual in- 
crease of the streamwise velocity with y. There is also a discrepancy in the 
data for the streamwise velocities at the reattachment point in Figure 4.2.2 and 
4.2.3. The cause of this discrepancy between the two data sets is not certain. 
However, comparison of the computed results with the data very close to the 
wall, y < 0.2 h, shows generally good agreement. 

Mean temperature profiles at various streamwise locations downstream of 
the step in Figure 4.2.4 show the steepest temperature gradients in the region 
very close to the wall. The large temperature gradient across the shear layer 
near the step far from the wall is generally well predicted by the present com- 
putations. The dark markers aty = 0 indicate the wall temperatures measured 
by Vogel and Eaton, and the vertical lines indicate numerical results of wall 
temperatures by the two methods, for x* = — 0.35, 0.05 and 0.45. The pre- 
diction by the low-Re model shows closer agreement to the experimental data 
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than that by the wall function method of Collins. Figure 4.2.5 shows an en- 
larged plot of the mean temperature profiles at x' = - 0.95 and —0.75. The 
wall temperature prediction by the low-Re model is again better than that of 
wall function method, except the large overprediction by the methods very 
near the step, at x* = — 0.95. 

Figure 4.2.6 shows the static pressure profiles on the top and bottom walls 
for Re e = 28,000 and 3je = 1.1 . The gradual increase of the static pressure 
on the top wall and the rapid increase in the bottom wall static pressure 
through reattachment to the downstream arc well predicted by the present 
computations. The computation by the low-Re model shows correctly the 
slightly accelerating flow in the upstream near the step. Pressure recovery 
downstream of reattachmcnt is predicted better by the low-Re model than by 
the wall function method. 

Figure 4.2.7 shows that the computed St at the reattachment point is about 
22% overpredicted by the low-Re model whereas the wall function method of 
Collins gives a 13% underprediction. The temperature difference, T w - 7^, 
at the reattachment point was about 4°C in the data of Vogel and Eaton, and 
if the reference condition is taken at 20°C , the maximum St is 0.00494 with 
their condition of U ref = 1 1.3 m/s and q w = 270 Wjm 2 . There is also a dis- 
crepancy in the profile of St for <5/e = 1.1 and Re e = 28,000 in their report. 
These uncertainties in the experimental data will have to be resolved. The 
present computation by the wall function method gives a peak St that is 13% 
lower than that reported by Collins, and the cause of this disagreement is not 
certain. The wall shear stress profile in Figure 4.2.7 shows that the computa- 
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tion by the low-Re model gives poor prediction in the recirculation region 
whereas prediction in the redeveloping region is much better than the wall 
function method. 
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4.3 FLOW IN CHANNEL WITH RECTANGULAR RIBS 


Parametric studies for repeated rib roughness were performed for 
ejD = 0.056 
e\b = 0.67, 1 
L\e =5, 7.5, 10, 15, 20 
Re = 5,200 — 41,800 

Air was chosen as the fluid, and properties were evaluated at 293 K. The wall 
heat flux used was 700 W j m 2 . A typical grid used for Lfe = 10 was NX = 84 
and NY = 74 with 60% of the grids located between the wall and the rib top. 

Figure 4.3.1 shows the velocity distribution for L\e = 5, 10 , and 20 at 
Re — 20,900. For Lje = 5, a large recirculation is present between the ribs and 
there is no flow reattachment. Though not shown in the figures, a small 
counterrotating vortex was observed in the corner behind the first rib. A large 
vortex near the second rib is shown in Figure 4.3. Id, and the velocity plot 
shows rapid acceleration of fluid along the front face of the rib. The static 
pressure contours for Lfe = 5 and 15 are shown in Figure 4.3.2. The turbulent 
kinetic energy contours in Figure 4.3.3 show maxima at the sharp edge of the 
second rib where flow impinges and a highly turbulent shear layer is generated. 
A second maximum occurs near the flow reattachment for Lje = 10, as ob- 
served in experiments. 

Isotherms for Lje = 5 and 10 in Figure 4.3.4 show high temperature regions 
where recirculations occur. Figure 4.3.5 shows the prediction of temperature 
along the channel wall including the front and back side of the rib for Lje = 
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5. The region between the first two vertical markers in the figures corresponds 
to the front side of the rib, and that between the second and third markers 
corresponds to the top of the rib. The temperature is the lowest at the top of 
the front side of the rib where there is a flow impingement, and the highest 
heat transfer coefficient is obtained. The boundary-layer development at the 
top of the second rib also leads to low wall temperature, and there is a sudden 
increase in temperature at the back side of the rib. Another maximum in local 
heat transfer coefficient between the ribs lies close to the second rib where the 
faster moving fluid above the ribs has penetrated. The temperature and local 
St distribution along the wall for Lje = 10 are shown in Figures 4. 3.6-9, and 
those for Lje = 15 are shown in Figures 4.3.10-13. The figures show the 
lowest wall temperature and the highest heat transfer coefficient between the 
ribs at the reattachment point. 

Figure 4.3.14 shows a comparison of the average friction factor with the 
experimental data of Han et al. [27]. The agreement of the friction factor with 
experiment is generally good, except for the underprediction of about 10% for 
both L/e = 10 and 15 at Re of about 20,000. For Lje = 5 and 7.5, well con- 
verged results for Re < 20,000 were difficult to obtain, and they are not pre- 
sented at this point. The prediction of the average St shows fair agreement 
with experimental data, except for generally low predictions. As shown in 
Figure 4.3.1, a large portion of area near the wall contains velocity vectors that 
are at some angle with the grid lines. Though the angles may be small, nu- 
merical diffusion effects are expected to be appreciable at higher Reynolds 
numbers. The heat transfer prediction is also sensitive to the details of near- 
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wall flow field such as the location and length of reattachment. Higher order 
models, such as the algebraic stress model, might give a superior simulation of 
this complex flow. A higher order differencing scheme such as QUICK 
method, may reduce the numerical error. 

Prediction of the average friction factor and Stanton number with each 
outer iteration are shown in Figure 4.3.15 for three cases. The average values 
for both /and St are stabilized after about 10 outer iterations. Figure 4.3.16 
is a typical plot of local Stanton number along the tube wall after each outer 
iteration. The total number of outer iterations shown is 20, and gradual con- 
vergence is obtained, though the behavior is erratic. 

Figure 4.3.17 shows a comparison of local Nusselt number distribution with 
the experimental data of Liou and Hwang [98]. The experiment was per- 
formed in a rectangular channel with an aspect ratio of 4:1, but the present 
computations were performed for a 2-dimensional channel. N in the Figure 
indicates the index of the rib in the heated section of the channel in their ex- 
periment. The comparison shows a fairly good prediction by the present re- 
sults except the front and top of the rib where large overprediction occurs. 


4.4 FLOW IN TUBE WITH RECTANGLAR RIBS 

Computation of flow in a channel with rectangular cross section is more 
appropriate in the design of a heat exchanger for active cooling, and thus the 
focus of the present work was on a 2-dimensional channel. The well known 
experiments by Webb et al. [26] using rectangular rib roughness were per- 
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formed in a tube, with parameters which were very different to those in the 
experiment of Han et al. [27], The tube had a larger pitch to rib height ratio 
and a smaller rib height to the channel height ratio compared to the channel 
of Han et al., and thus presented a different test for the turbulence model. 
This problem was actually investigated first, and is included here for com- 
pleteness. 

Parametric studies for repeated rib roughness were performed for 
e/D = 0.02 
ejb = 1 

Lie =20 and 40 
Re= 18,940 - 56,810 

The ratio of roughness height to width, e\ b , was 1.94 in the experiments of 
Webb et al. [26]. For the large pitch to roughness height ratio, this variation 
of e/b should have negligible effect on pressure drop and overall heat transfer 
rate in the fully rough regime. With the outlet at the edge of a rib, the wider 
rib shape of ejb = 1.0 gave more stable numerical solutions. Air was chosen 
as the fluid, and properties were evaluated at 300 K. The wall heat flux used 
was lOOWjm 2 . The typical grid used for L/e = 20 was NY = 85, NX = 84 
with 30 grids in the radial direction between the wall and rib tip. 

Figure 4.4.1 show the prediction of temperature and the local Stanton 
number distribution along the pipe wall for Lje = 20. Results are shown be- 
tween two slabs before the first rib and the slab NX - 1 on the second rib. 
Near-symmetry of wall temperature and Stanton number was obtained be- 
tween the inlet and outlet. The maximum temperature occured in the corner 


81 



of the first rib. The maximum in local Stanton number between two ribs oc- 
curs at the flow reattachmcnt point. A second maximum heat transfer occurs 
at the sharp edge of the second rib where there is flow impingement and 
boundary-layer development. 

Figure 4.4.2 show a comparision of the average friction factor and Stanton 
number with experimental data of Webb ct.al. [26]. The agreement with ex- 
periment was generally quite good, except for the heat transfer prediction for 
L I e = 40. A stabilized heat transfer coefficient was difficult to obtain at high 
Reynolds number for this pitch over roughness height ratio. 


82 



4.5 JET IMPINGEMENT COOLING 

Numerical computations were performed for turbulent axisymmetric jets 
impinging on a heated flat plate for various values of Re and H/d. The 
Reynolds number was varied from 23,750 to 80,400, and H/d was varied from 
2 to 22. The typical grid of the computation domain for Hid = 6 was 58X74 
cells in the x and r directions with about 44% of the grids located within the 
normal distance of one jet diameter from the impingement wall and dense grids 
located around the stagnation point. The outlet was located at either 10 or ,20 
jet diameters downstream from the stagnation point depending on the distance 
from the jet outlet to the impingement surface. Either fully developed profiles 
for velocity, turbulent energy, and dissipation rates obtained from the numer- 
ical solution of the turbulent flow in a pipe, or uniform profiles were used at 
the jet outlet. 

Figure 4.5.1 shows the velocity vector plots for Hjd = 2 and 6 computed 
with a low-Re turbulence model and Yap's correction. The jet grows by 
entrainment of the stagnant fluid and deflects into a radial wall jet. The decay 
of the jet centerline velocity in Figure 4.5.2a shows that the jet flow is not 
much affected by the impingement wall up to about one jet diameter from the 
wall. About 95% of the flow deceleration occurs within one nozzle diameter 
from the wall for Hjd = 6. Comparison of the jet centerline velocity decay 
with the experimental data of Giralt et al. [43] for several nozzle heights is 
shown in Figure 4.5.2b. Rapid decay of the centerline velocity to zero at the 
stagnation point is shown for Hjd = 6.67. Computations for larger nozzle 
heights show characteristics of the free jet data, and predicts correctly the in- 
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fluence of the impingement wall on the oncoming Ao\v(which extends further 
upstream as the nozzle height increases). The slightly faster decay near the 
impingement wall for Hjd = 22 as compared to the experimental data could 
be due to greater spreading rate of the free jet. 

Figure 4.5.3a shows the computed results for the axial velocities, uju m plot- 
ted against the nondimensional parameter, rj5 u , where 5 U is the radial dis- 
tance r where u = 1/2 u m . Similar proAles at different axial locations compare 
fairly well with the experimental data, except farther away from the center of 
the jet, for rj5 u > 1.5. The experimental data of Beltaos and Rajaratnam 
[42] were found to agree well with the theoretical solution of Tollmien [70] 
that is based on Prandtl's mixing-length hypothesis. Tollmien's solution is 
known to describe accurately the velocity distribution in a two-dimensional 
turbulent jet, and thus the present numerical computation seems to give fairly 
good prediction of Aow Field in the jet development region. 

The distribution of the axial velocities at various locations in Figure 4.5.3b 
shows the continuous decay of the centerline velocity as the impingement wall 
is approached. Though not plotted in the figure, the radial velocity near the 
center of the jet is positive causing the jet to spread out, and the radial velocity 
further from the center is negative and Aow entrainment occurs. 

Dimensionless profiles of mean radial velocities are shown in Figure 4.5.4. 
V m is the maximum velocity at each station and 5 V is the height where 
V = 0.5 V m . Comparison with experimental data of the classical wall jet ob- 
tained by Schwarz and Cosart [99] is in good agreement except that the lo- 
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cation of the maximum radial velocities arc slightly underpredicted. Only 
portion of the experimental data are plotted in the figure. 

The distribution of the static pressure along the impingement wall for Hjd 
= 6 is shown in Figure 4.5.5. Due to the uncertainty of the outlet condition 
of the jet in the experiment by Giralt et al. , two different jet outlet profiles, one 
which was obtained numerically at 50 diameters downstream in a tube, and 
one with uniform profiles, were employed. The uniform jet nozzle outlet con- 
dition fits the experimental data more closely than the fully developed condi- 
tion. Giralt et al. investigated the pressure distributions over the range 1.2 
< Hjd < 20 at 30,000 < Re < 80,000 and narrower distribution of the wall 
static pressure was obtained with increasing jet height, for Hjd < 8. They at- 
tributed this to the more concentrated kinetic energy of the jet near the 
centerline with more non-uniform velocity profiles at the start of the jet 
impingement region. The present computations with a fixed jet height and two 
different jet outlet conditions show the same trend. About 35% increase in the 
wall shear stress due to probably the steeper pressure gradient for the fully 
developed jet nozzle condition is shown in Figure 4.5.5b. This trend is con- 
sistent with Amano and Brandt [50] whose prediction for Hjd = 18 and Re 
= 180,000 show 26% increase in wall shear stress when 5.5 power velocity 
profile was used. Amano attributed this to the higher turbulent kinetic energy 
at the edgy of the jet in the case of 5.5 power velocity profile. Figure 4.7.5b 
also shows a shift in the location of the maximum shear stress toward the 
stagnation point, a decrease of about 18%. 
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Comparison with the experimental data of Baltaos and Rajaratnam is made 
for Hid = 21.2 and Re = 80,400 in Figure 4.5.6a. The pressure distribution 
along the impingement wall show a fairly good agreement except for overall 
ovcrprediction. The radial location where the pressure is 10% ot the stag- 
nation pressure is at about 3.6 d compared to about 1.1 d for Hjd = 6 as 
shown in Figure 31. a. The centerline static pressure in Figure 4.5.6b shows 
good agreement with the data of Beltaos. Predictions for both H/d = 6 and 
21.2 show no effect of the impingement wall on the jet up to 0.8 H down- 
stream. For Hjd = 21.2, computation shows slightly negative static pressure 
indicating that pressure in the free-jet region is below the ambient pressure. 

The wall shear stress prediction in Figure 4.5.7a shows about 31% over- 
prediction in the maximum shear stress, whereas good agreement is achieved 
away from the stagnation region, r/// > 0.15. Comparison with more exper- 
imental data is needed to resolve the discrepancy in the stagnation region. 
Baltaos and Rajaratnam indicate that the error in the shear stress for r/H less 
than about 0.08 could be greater than 6% due to the presence of large pressure 
gradients. Prediction of the maximum radial velocity along the impingement 
wall in Figure 4.5.7b shows that the peak is only about 22% of the jet outlet 
velocity. 

Figure 4.5.8a shows the half-width spreading rate of the round free jet for 
Hjd = 21.2 and Re = 74,000. Numerical results for 0.4 < x/H < 0.9 predicts 
the spreading rate, dS u / dx, equal to 0.103, compared to 0.093 of the data of 
Baltaos. Malin [54] reported the spreading rate of round free jet equal to 
0.113 with the standard k-e model. The half-width spreading rate of the 
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radial wall jet is shown in Figure 4.5.8b. The spreading rate, dd Y j dr, for 
Hid = 10 and Re = 23750 was 0.09, and that for Hid = 21.2 and Re =74000 
was 0.085 over the range of 0.4 < r/// < 1.0 . The combined correlation for 
both Hjd = 10 and 21.2 indicate the average spreading rate of 0.09, or <5 V / H 
is expressed as 

b v l H = 0.0948 (r / //) 0 - 903 (167) 

Porch et al. [40] investigated turbulent radial wall jet for 8 < H / d < 24 and 
64,000 < Re < 288,000 ,and the increase of jet thickness for 0.5 < r/H< 3 
was correlated by 

S v l H = 0.098 (r/ H) 09 (168) 

Comparing the two correlations shows that the prediction by the low-Re model 
is fairly good, and the average spreading rate of 0.9 compared to the exper- 
imentally excepted value of 0.085 - 0.095 is superior to the computed value of 
0.068 by the standard k — c model. 

Similarity profiles for k and uv in a free jet computed for Hjd — 21.2 show 
fair agreement with the experimental data of Wygnanski and Fiedler [44] in 
Figure 4.5.9. The uv profile has the maximum value at about r = 0.65 b u 
from the center of the jet where the velocity gradient dU I dr is the largest. The 
present results are also in very close agreement with the numerical results of 
the standard k — e model that was reported by Malin who computed free and 
wall jets with various modifications to k — e and k — IV model. Figure 4.5.10 
shows computed profiles of the turbulent kinetic energy and the turbulent 


87 



shear stress for H I d = 10 and 21.2 at various radial locations. Computed re- 
sults of the turbulent intensity is underpredicted by about 17% compared to 
the data of Ng [100], but Malin points out that the data of Ng cannot be re- 
garded as definitive, and that Poreh ct al. measured a peak intensity of 
£ 1 / 2 1 = o.35. (The experimental data of Ng was reproduced from the pa- 

per of Malin due to unavailability of Ng's thesis at present). The predicted 
profiles for H j d = 10 is very close to the numerical results of the standard 
k - e model that was reported by Malin, except that the present computat'pn 
of uv with the low-Re model gives up to 20% improvement near the wall, 
x < 0.5 <5 V . in the case of the radial wall jet, the similarity profiles for k and 
uv are not quite obtained for x > 0.5 (5 V , and since only one computed profile 
for each case is given from Malin, it was not possible to compare the degree 
of similarity. However, the experimental data of Poreh et al. show scatter in 
measurement of both k and uv, and they summarized that it was difficult to 
conclude from the data whether the turbulent shear stress in the wall jet is 
similar. 

Figure 4.5.1 1 shows the Nusselt number profile for H / d = 2, 6, 10, 14, and 
Re = 23,750. The outlet profiles for the jet were obtained from the numerical 
solution from a pipe with the length of 72 diameters to simulate the exper- 
imental condition. Comparison with the experimental data of Baughn and 
Shimizu shows poor agreement near the stagnation point, as much as 56% 
overprediction for Hj d = 6. In the case of Hj d = 2, the minimum Nu at 
rfd of about 2 are well produced by the computation. The predicted profiles 
are in good agreement with the data for the downstream region. Figure 4.5.13 
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and 4.5.14 show the results when the extra source term in the e equation is 
damped more than in the origional proposal of Yap. The predicted profiles 
of the local Nu show significant improvements except for the absence of the 
second peak for Hid = 2. The extra source term had to be about 1.6 to 6 
times larger depending on H / d. More detailed investigation is needed to 
produce the necessary correction term to improve the heat transfer prediction 
in the vicinity of the stagnation point. 
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4.6 VARIABLE PROPERTIES CHANNEL FLOW USING THE 
LOW-RE MODEL 

The effect of variable properties on heat transfer augmented by ribs in 
channel flow can be deduced from empirical equations, such as that of Vilemas 
and Simonis [38], 

Nu b = 0.029 Re b M Pr b 6 (-^-) ( 1 69) 

1 b 

where n depends on geometry and Reynolds number, and is given by Eq.(29). 
Thus 

Sf 4 «4' 16 (2?4” ; n<0 (170) 

J b 

With heating n b increases and TJT b decreases in the axial direction, thus both 
terms contribute to the increase of St b along the channel. 

The friction factor for flow over sand-grain roughness is known to be func- 
tion of geometry only, and independent of Reynolds number in the fully rough 
regime where form drag is the dominant mechanism. But, the friction factor 
for flow over repeated-rib roughness involves more geometric parameters, 
namely the height of roughness e, the width of the roughness w, and the pitch 
of roughness elements L. For widely spaced ribs, the flow reattaches behind 
each rib, and a viscous layer grows. This viscosity dependent shear stress 
should be dominant as L\e approaches infinity. Savage and Myers [101] 
performed experiments with water flow in a tube with rectangular rib 
roughness. They investigated the contribution of form drag to the total forces 
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retarding flow for different roughness configurations and Reynolds numbers. 
For Lje = 1, the total pressure drop was entirely due to a skin-friction effect. 
The relative contribution of skin friction decreased almost linearly as Lje in- 
creased and appeared to pass through a minimum at about Lje = 13. For 
ejD = 0.04 and Re= 10 5 , the minimal contribution of skin friction was 30 % 
of the total retarding force at Lje ^ 10. For the repeated-rib roughness, it 
seems that form drag may not be the dominant mechanism and skin friction 
may not be negligible, depending on the spacing of ribs. Consequently »,he 
friction factor may depend on viscosity in a variable properties flow even 
though the equivalent sandgrain roughness indicates the fully rough regime. 
The variation of friction factor downstream along the channel wall may be 
difficult to correlate since the relative contribution of skin friction depends on 
the rib spacing. 

As a preliminary study, a channel with Lje = 15 , e/Z) = 0.056 was chosen 
to observe the effect of variable physical properties. The low- Re model with 
the Yap's modification was employed. The inlet Reynolds number was varied 
from 11,500 to 20,000, and the calculated temperature ratio, TJT b , was be- 
tween 1.42 and 1.88. Pr and c p were assumed constant for these initial calcu- 
lations. Iterations were performed until a stabilized value of friction factor and 
Stanton number were obtained. In the present computation, the calculated 
variables at the outlet after an outer iteration were substituted as inlet condi- 
tions for the next outer iteration. Thus, the temperature level is increasing in 
the computation domain for successive outer iterations, and the process is very 
similar to marching in the flow direction. A slight increase in / and St were 
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observed in the course of iteration at high heat fluxes, and it is suspected that 
the heat transfer coefficient may not have stabilized due to flow laminarization 
in this low Reynolds number range. 

Figure 4.6.1 shows a plot of fjf c versus TjT b for Re = 10,000 and 15,000, 
where / c is the constant property value. The Reynolds number is an average 
value calculated with bulk properties over a pitch. The present calculation 
shows about 13% and 21% reduction in /at Re = 15,000 and 10,000 respec- 
tively, when the air was heated to about T w jT b = 1.7. The results also shQW 
an increase in / with increasing Re , and the effect of Re on /seems not negli- 
gible. The experimental data of Vilemas and Simonis [38] showed a large ef- 
fect of temperature ratio on friction at low range of Re but smaller effect of 
Re on friction. They found a 20% decrease in /in the Reynolds number range 
of 10,000 to 20,000 when the channel was heated to TJT b = 1.8 from an 
adiabatic flow condition. The effect of Reynolds number on /was negligible 
for heating up to TJT b = 1.8 in the range of Reynolds number less than about 
60,000. However the data for the influence of temperature ratio on friction is 
for one channel only (e/Z) = 0.013, L/e = 9.5), and the effect of the parameter 
eld e could not be discerned. The friction factor measured in the experiment 
of Vilemas and Simonis was the overall friction factor with an outer smooth 
wall and inner rough wall, whereas the present computation is performed with 
the both walls of identical roughness. Thus an exact comparison cannot be 
made. The channel of the present computation also has much larger pitch to 
roughness height ratio than the test section of Vilemas and Simonis, and the 
ratio of the roughness height to the channel height of flow passage is about 
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11% compared to only 2.6% in the experiment. Thus, the channel of the 
present study has larger flow recirculation and redeveloped region than the 
experimental setup. 

Figure 4.6.2 shows the plot of StjSt c versus T w j T b . There is about 5% re- 
duction in heat transfer rate for both Re = 10,000 and 15,000 when the air is 
heated to T w / T b ~ 1 .7. There is no apparent effect of Reynolds number on 
heat transfer for the small range of Reynolds number in the present computa- 
tion. This seems to be supported by the experimental results of Vilemas a*?d 
Simonis which show only about 5% difference in Nu / Nu c from Re = 1 1,000 
to 31,000 when the air was heated to T w j T b ~ 1.7. However, their results 
show a large reduction in heat transfer with heating in the low Reynolds 
number range: about 19% at Re = 11,000 when the gas was heated to 

T w I T b ~ 1 .7. More computations need to be performed in order to investigate 
the effect of Reynolds number and the parameter ejD e . 


4.7 VARIABLE PROPERTIES FLOW WITH WALL FUNCTIONS 

In the present computations, Pr and c p are assumed constant, and p and p 
are assumed to vary as 



0.67 


= (^-) 


in 



(171) 


The wall function approach for variable properties flow described in section 
2.2.5 was used to calculate the flow in channel with ribs, and the prediction 
of heat transfer showed opposite trend with heating: an increase of heat 
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transfer with temperature. As stated by Viegas and Rubesin, it seems that 
improvements in eddy viscosity modeling and careful consideration of the 
mean values of velocity and temperature in the first mesh volume arc required 
in case of large wall heat fluxes. But, more importantly it seems that the in- 
crease of the viscous sublayer with heating should be incorporated into the 
Ciofalo and Collins approach. The proposed approach is to include the tem- 
perature ratio in case of variable properties flow as below. 



where i f/ e and y vo are the turbulence intensity and the non-dimensional sublayer 
thickness respectively in the equilibrium boundary layer, and a and b are em- 
pirical constants. This modification will make the P - function in the law of 
the wall for temperature a function of heating. Even though, this is possible 
with experimental data of variable properties flow, it would be similar to the 
derivation of the standard wall functions. Another drawback of the wall 
function approach was an inability to predict the recirculation zone at the top 
of the forward-facing rib, causing a low pressure zone above the rib which in 
turn gave a larger overall pressure drop. 

Though this is less important, a possible improvement over Eq. (97) would 
be to adopt Amano's approach of evaluating the production and destruction 
terms in the £ equation, taking into account the variation of k in the near wall 
layer. A simpler approach would be to use Eq. (95) as a mean value. These 
possible improvements need to be further investigated. 
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4.8 EFFECT OF WALL THERMAL BOUNDARY CONDITIONS. 

The computations presented were all performed with a constant wall heat 
flux boundary condition to correspond with the majority of the experiments 
that were used for comparison purposes. Subsequently flow in a channel with 
transverse rectangular ribs was solved with a constant wall temperature 
boundary condition, and an approximately four-fold increase of heat transfer 
coefficient at the top edge of the forward-facing rib was found. In the recir- 
culating region behind the rib, the prediction of heat transfer was comparable 
with that for a constant wall heat flux boundary condition. 

In order to investigate this unexpected results further, a simpler problem 
was chosen, namely a laminar flow over a forward-facing step. Figure 4.8.1 
shows a velocity field for Re = 550 and e j H = 0.13. The grid employed was 
52x28 with 10 grid nodes placed below the step. Figure 4.8.2 shows the local 
Nusselt number distributions along the channel wall including the front face 
of the step. For the both boundary conditions, the Nusselt number reaches a 
minimum in the corner of the step, and there is a gradual increase in Nusselt 
number along the front face of the step reaching a maximum on the top of the 
step. For a constant T w , the minimum is much smaller and the maximum 
much larger(2 times) than for a constant q w boundary condition. 

Grid refinement tests are shown in Figure 4.8.3 for four different grid ar- 
rangements. The solutions are practically identical for all cases except for the 
overprediction of heat transfer coefficient on the top of the step with the grid 
of 40x20. The prediction over the face of the step is invariant for all the grid 
arrangements. 


95 



Isotherms for the both wall boundary conditions are plotted in Figure 4.8.4 
where the scale of y coordinate has been enlarged to have a better view of 
temperature profiles below the step. The temperature contour has been plotted 
from the inlet to the forward-facing step, and from the wall to the top of the 
step only. The solution with the constant wall temperature boundary condi- 
tion produces a steeper temperature gradients along the front face of the rib 
compared to that with wall heat flux boundary condition. In the corner of the 
step, the velocity of fluid is very small. Thus the fluid next to the wall retains 
a temperature very close to the w T all temperature in the case of constant T w 
boundary condition, and the fluid acts as an insulating layer. This leads to the 
very low heat transfer coefficient in the corner of the step as shown in Figure 
4.8.2. The constant q w boundary condition establishes wall temperature gra- 
dient along the face of the step, and this leads to steeper temperature gradient 
of fluid layer along the step wall , and subsequently higher heat transfer coef- 
ficient in the corner of the step than that for the constant T w boundary condi- 
tion. Whereas the wall temperature at the top of the step for the constant q w 
boundary condition has cooled to about 80% of the wall temperature at the 
base of the step, the wall temperature along the front face of the step for the 
constant T w boundary condition is fixed at a constant value. This effect seems 
to be similar to the entrance region problem where the temperature gradient 
at the wall is theoretically infinite. In fact, the temperature gradient at the top 
edge of the step is almost identical to that at the wall of the second slab from 
the entrance in the present computation domain. These effect of wall bound- 
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ary conditions on heat transfer was not observed in cither abrupt pipe expan- 
sion flow or in jet impingement cooling. 
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Chapter V 
CONCLUSIONS 

5.1 OVERALL CONCLUSIONS 

The flow downstream of a sudden pipe expansion and over a backward- 
facing step were solved with both high and low Reynolds number k - e model. 
The standard wall function method proved to be inadequate in the calculation 
of heat transfer rates for flow with reattachmcnt. The wall function approach 
of Ciofalo and Collins that relates the viscous sublayer thickness to the near- 
wall turbulence intensity gave better heat transfer prediction than the standard 
approach, but the overall performance was poor. The best prediction of heat 
transfer was obtained with the Yap modified Jones and Launder low Reynolds 
number k — e model. 

The large part of computations were performed on two elliptic flows: the 
flow in a channel with rectangular transverse ribs and jet impingement on a 
flat plate. The objectives were the application of turbulence models to these 
flows to facilitate the design of active cooling in hypersonic flight, and further 
improvement in the turbulence model for better prediction of heat transfer. 

A simple periodicity boundary condition was used to compute the fully de- 
veloped flow in a channel with repeated rib roughness. The low- Reynolds 
number k — e model with Yap's modification gave an adequate agreement with 
experiment for constant property flow. The results for average heat transfer 
rates were generally lower than the experimental data. For jet impingement 
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cooling, there were overprediction of heat transfer rates in the neighborhood 
of reattachment, and modification to the extra source term in the dissipation 
rate equation was required. 

Investigation of the variable property flow in a channel with ribs with the 
same turbulence model showed reduction in both friction and average heat 
transfer rates with heating. The effect of Reynolds number on the friction was 
not negligible in contrast to heat transfer rates. The computation domain of 
the present study was different from the test section of the experiment, and 
exact conclusion could not be made. The wall function method was modified 
to account for the effect of variable properties with non-adiabatic surface, and 
results showed that a further near-wall modelling is required to incorporate the 
increase of the viscous sublayer with heating. 

There was a large increase in heat transfer rate at the top edge of the rib 
when the wall boundary condition was changed to constant wall temperature 
from constant wall heat flux. A laminar flow over a forward-facing step was 
chosen to further investigate the effect of wall thermal boundary condition. 
The Nusselt number at the edge of the step for the constant T w boundary 
condition was about two times higher than that for the constant q w boundary 
condition. 
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5.2 RECOMMENDATIONS 

Various turbulence models have been shown to perform fairly well in pre- 
diction of mean flow by many investigators, but there needs to be a more in- 
tensive investigation to improve the turbulence models for better prediction of 
heat transfer in complex flows. More experimental measurements of turbu- 
lence quantities and heat transfer rates are needed especially for the stagnation 
point region of the jet impingement flow, and for variable property flow with 
heating. These will facilitate for better modelling of turbulence models. 
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Figure 3.2 The staggered control volume in the x-v plane. 
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Figure 4.1.1 Velocity profiles for (a) djD = 0.4, Re = 12,310 (b) d\D — 0.8, 
Re = 20,130. 
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Figure 4.1.2 Comparison of Nusselt number and skin-friction coefficient with 
the data of Baughn et al. djD = 0.4, Re = 12,310. 
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Figure 4.1.3 Comparison of Nusselt number for (a) Re = 23,210, (b) Re = 
40,750 with the data of Baughn et al. 
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Figure 4. 1 .4 Comparison of Nusselt number and skin-friction coefficient with 
the data of Baughn ct al. dj D = 0.8, Re = 10,570. 
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Figure 4.1.5 Comparison of Nusselt number for djD = 0.8 (a) Re = 20,130 
(b) Re =39,300. 




Figure 4. 1 .6 Comparison of Nusselt number and skin-friction coefficient with 
the data of Baughn et al. for d\D = 0.8, Re = 76,080. 




Figure 4.2.1 Vector plot for flow over a backward-facing step, expansion 
ratio =1.25. Re e — 28,000. 
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Figure 4.2.3 Mean velocity profile in the near-wall region. Re e = 28,000. 



Figure 4.2.4 Temperature profile for flow over a backward-facing step. Re e 
= 28,000. 



Figure 4.2.5 Temperature profile in the near-wall region. Re e = 28,000. 
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Figure 4.2.6 Wall static pressure profile for flow over backward-facing step. 
Re e = 28,000. 



Figure 4.2.7 Predictions of Stanton 
backward-facing step. Re e = 28,000. 





















Figure 4.3.3 Turbulent kinetic energy contour for Re = 20,900. (a) Lje = 
5, (b) Lje = = 10. 
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Figure 4.3.9 Prediction of wall temperature and Stanton number for Lje 
10, Re = 41,800. 
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Figure 4.3. 
15, Re = 2 
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Figure 4.3.13 
15, Re = 41,8 
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Figure 4.3.15 Average friction factor and Stanton number with outer iter- 
ation. 
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Figure 4.5.1 Vector plots for axisymmetric jet impingiment on a flat surface 
for Re = 23,750. (a) Hjd = 2, (b) Hjd = 6. 
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Figure 4.5.2 Comparison of the centerline velocity decay for (a) H\d = 2 and 
6 at Re = 23,750. (b) H/D = 6.67 and Re = 23,750. H)D = 15.56 and 
Re = 40,000. HID = 22 and Re = 80,400. 








Figure 4.5.4 Radial wall jet profile for Hjd = 6 and Re = 23,750. 
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Figure 4.5.5 Comparison for Hjd = 6 and Re = 30,000. 
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Figure 4.5.6 Comparison for Hid = 21.2 and Re = 80,400. (a) wall static 

pressure (b) centerline static pressure. 





Figure 4.5.7 Predictions for Hjd = 21.2 and Re = 80.400. (a) wall shear 

stress (b) maximum radial velocity. 
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Figure 4.5.8 Predictions of spreading rates for Hid =21.2 and Re = 74,000, 
Hjd = 10 and Re = 23,750. (a) free jet (b) wall jet 
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Figure 4.5.9 Prediction of (a) free-jet k profile, (b) free-jet uv profile for Re 
= 74,000. 






Figure 4.5.10 Prediction for Hjd =21.2 and Re = 74,000, Hjd = 10, and 
Re = 23,750. (a) wall-jet k profile, (b) wall-jet uv profile 
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Figure 4.5.12 Comparison of Nussclt number profile with modified damping 
in £ equation. 



Figure 4.5.13 Comparison of Nussclt number profile with modified damping 
in e equation. 
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